Belle II Software  release-05-01-25
EnergyLossForExtrapolator.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010-2011 Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Leo Piilonen *
7  * Derived from: G4EnergyLossForExtrapolator.cc *
8  * (use geant4e-specific particles; include pions & kaons) *
9  * *
10  * This software is provided "as is" without any warranty. *
11  **************************************************************************/
12 
13 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
14 
15 #include <simulation/kernel/EnergyLossForExtrapolator.h>
16 
17 #include <CLHEP/Units/PhysicalConstants.h>
18 #include <CLHEP/Units/SystemOfUnits.h>
19 
20 #include <G4PhysicsLogVector.hh>
21 #include <G4ParticleDefinition.hh>
22 #include <G4Material.hh>
23 #include <G4MaterialCutsCouple.hh>
24 #include <G4ParticleTable.hh>
25 #include <G4LossTableBuilder.hh>
26 #include <G4MollerBhabhaModel.hh>
27 #include <G4BetheBlochModel.hh>
28 #include <G4Version.hh>
29 #if G4VERSION_NUMBER < 1001
30 #include <G4eBremsstrahlungModel.hh>
31 #else
32 #include <G4eBremsstrahlungRelModel.hh>
33 #endif
34 #include <G4MuPairProductionModel.hh>
35 #include <G4MuBremsstrahlungModel.hh>
36 #include <G4ProductionCuts.hh>
37 #include <G4LossTableManager.hh>
38 #include <G4WentzelVIModel.hh>
39 
40 #include <framework/logging/Logger.h>
41 
42 using namespace Belle2;
43 using namespace Belle2::Simulation;
44 
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
46 
47 EnergyLossForExtrapolator::EnergyLossForExtrapolator(void) :
48  m_Particle(nullptr),
49  m_Electron(nullptr),
50  m_Positron(nullptr),
51  m_MuonPlus(nullptr),
52  m_MuonMinus(nullptr),
53  m_PionPlus(nullptr),
54  m_PionMinus(nullptr),
55  m_KaonPlus(nullptr),
56  m_KaonMinus(nullptr),
57  m_Proton(nullptr),
58  m_AntiProton(nullptr),
59  m_Deuteron(nullptr),
60  m_AntiDeuteron(nullptr),
61  m_ProductionCuts(nullptr),
62  m_DedxElectron(nullptr),
63  m_DedxPositron(nullptr),
64  m_DedxMuon(nullptr),
65  m_DedxPion(nullptr),
66  m_DedxKaon(nullptr),
67  m_DedxProton(nullptr),
68  m_DedxDeuteron(nullptr),
69  m_RangeElectron(nullptr),
70  m_RangePositron(nullptr),
71  m_RangeMuon(nullptr),
72  m_RangePion(nullptr),
73  m_RangeKaon(nullptr),
74  m_RangeProton(nullptr),
75  m_RangeDeuteron(nullptr),
76  m_InvRangeElectron(nullptr),
77  m_InvRangePositron(nullptr),
78  m_InvRangeMuon(nullptr),
79  m_InvRangePion(nullptr),
80  m_InvRangeKaon(nullptr),
81  m_InvRangeProton(nullptr),
82  m_InvRangeDeuteron(nullptr),
83  m_MscatElectron(nullptr),
84  m_Material(nullptr),
85  m_MaterialIndex(0),
86  m_ElectronDensity(0),
87  m_RadLength(0),
88  m_Mass(0),
89  m_ChargeSq(0),
90  m_KineticEnergy(0),
91  m_Gamma(1.0),
92  m_BetaGammaSq(0),
93  m_BetaSq(0),
94  m_Tmax(0),
95  m_LinLossLimit(0.01),
96  m_UserTmin(1.*CLHEP::MeV),
97  m_UserTmax(10.*CLHEP::TeV),
98  m_UserMaxEnergyTransfer(DBL_MAX),
99  m_Nbins(70),
100  m_NMaterials(0),
101  m_Initialised(false)
102 {
103  Initialisation();
104 }
105 
106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
107 
108 EnergyLossForExtrapolator:: ~EnergyLossForExtrapolator()
109 {
110  m_DedxElectron->clearAndDestroy(); delete m_DedxElectron;
111  m_DedxPositron->clearAndDestroy(); delete m_DedxPositron;
112  m_DedxMuon->clearAndDestroy(); delete m_DedxMuon;
113  m_DedxPion->clearAndDestroy(); delete m_DedxPion;
114  m_DedxKaon->clearAndDestroy(); delete m_DedxKaon;
115  m_DedxProton->clearAndDestroy(); delete m_DedxProton;
116  m_DedxDeuteron->clearAndDestroy(); delete m_DedxDeuteron;
117  m_RangeElectron->clearAndDestroy(); delete m_RangeElectron;
118  m_RangePositron->clearAndDestroy(); delete m_RangePositron;
119  m_RangeMuon->clearAndDestroy(); delete m_RangeMuon;
120  m_RangePion->clearAndDestroy(); delete m_RangePion;
121  m_RangeKaon->clearAndDestroy(); delete m_RangeKaon;
122  m_RangeProton->clearAndDestroy(); delete m_RangeProton;
123  m_RangeDeuteron->clearAndDestroy(); delete m_RangeDeuteron;
124  m_InvRangeElectron->clearAndDestroy(); delete m_InvRangeElectron;
125  m_InvRangePositron->clearAndDestroy(); delete m_InvRangePositron;
126  m_InvRangeMuon->clearAndDestroy(); delete m_InvRangeMuon;
127  m_InvRangePion->clearAndDestroy(); delete m_InvRangePion;
128  m_InvRangeKaon->clearAndDestroy(); delete m_InvRangeKaon;
129  m_InvRangeProton->clearAndDestroy(); delete m_InvRangeProton;
130  m_InvRangeDeuteron->clearAndDestroy(); delete m_InvRangeDeuteron;
131  m_MscatElectron->clearAndDestroy(); delete m_MscatElectron;
132  delete m_ProductionCuts;
133  for (const G4MaterialCutsCouple* couple : m_Couples) delete couple;
134  m_Couples.clear();
135 }
136 
137 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
138 
139 G4double EnergyLossForExtrapolator::EnergyAfterStep(G4double kinEnergy,
140  G4double stepLength,
141  const G4Material* mat,
142  const G4ParticleDefinition* part)
143 {
144  //- if (!isInitialised) Initialisation();
145  G4double kinEnergyFinal = kinEnergy;
146  if (SetupKinematics(part, mat, kinEnergy)) {
147  G4double step = TrueStepLength(kinEnergy, stepLength, mat, part);
148  G4double r = ComputeRange(kinEnergy, part);
149  if (r <= step) {
150  kinEnergyFinal = 0.0;
151  } else if (step < m_LinLossLimit * r) {
152  kinEnergyFinal -= step * ComputeDEDX(kinEnergy, part);
153  } else {
154  G4double r1 = r - step;
155  kinEnergyFinal = ComputeEnergy(r1, part);
156  }
157  }
158  return kinEnergyFinal;
159 }
160 
161 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
162 
163 G4double EnergyLossForExtrapolator::EnergyBeforeStep(G4double kinEnergy,
164  G4double stepLength,
165  const G4Material* mat,
166  const G4ParticleDefinition* part)
167 {
168  //- if (!isInitialised) Initialisation();
169  G4double kinEnergyFinal = kinEnergy;
170 
171  if (SetupKinematics(part, mat, kinEnergy)) {
172  G4double step = TrueStepLength(kinEnergy, stepLength, mat, part);
173  G4double r = ComputeRange(kinEnergy, part);
174 
175  if (step < m_LinLossLimit * r) {
176  kinEnergyFinal += step * ComputeDEDX(kinEnergy, part);
177  } else {
178  G4double r1 = r + step;
179  kinEnergyFinal = ComputeEnergy(r1, part);
180  }
181  }
182  return kinEnergyFinal;
183 }
184 
185 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
186 
187 G4double EnergyLossForExtrapolator::TrueStepLength(G4double kinEnergy,
188  G4double stepLength,
189  const G4Material* mat,
190  const G4ParticleDefinition* part)
191 {
192  G4double res = stepLength;
193  //- if (!isInitialised) Initialisation();
194  if (SetupKinematics(part, mat, kinEnergy)) {
195  if (part == m_Electron || part == m_Positron) {
196  G4double x = stepLength * ComputeValue(kinEnergy, m_MscatElectron);
197  if (x < 0.2) res *= (1.0 + 0.5 * x + x * x / 3.0);
198  else if (x < 0.9999) res = -std::log(1.0 - x) * stepLength / x;
199  else res = ComputeRange(kinEnergy, part);
200 
201  } else {
202  res = ComputeTrueStep(mat, part, kinEnergy, stepLength);
203  }
204  }
205  return res;
206 }
207 
208 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
209 
210 G4bool EnergyLossForExtrapolator::SetupKinematics(const G4ParticleDefinition* part,
211  const G4Material* mat,
212  G4double kinEnergy)
213 {
214  if (!part || !mat || kinEnergy < CLHEP::keV) return false;
215  //- if (!isInitialised) Initialisation();
216  G4bool flag = false;
217  if (part != m_Particle) {
218  flag = true;
219  m_Particle = part;
220  m_Mass = part->GetPDGMass();
221  G4double q = part->GetPDGCharge() / CLHEP::eplus;
222  m_ChargeSq = q * q;
223  }
224  if (mat != m_Material) {
225  G4int i = mat->GetIndex();
226  if (i >= m_NMaterials) {
227  B2WARNING("EnergyLossForExtrapolator: index i= " << i << " is out of table - NO extrapolation");
228  } else {
229  flag = true;
230  m_Material = mat;
231  m_ElectronDensity = mat->GetElectronDensity();
232  m_RadLength = mat->GetRadlen();
233  m_MaterialIndex = i;
234  }
235  }
236  if (flag || kinEnergy != m_KineticEnergy) {
237  m_KineticEnergy = kinEnergy;
238  G4double tau = kinEnergy / m_Mass;
239 
240  m_Gamma = tau + 1.0;
241  m_BetaGammaSq = tau * (tau + 2.0);
242  m_BetaSq = m_BetaGammaSq / (m_Gamma * m_Gamma);
243  m_Tmax = kinEnergy;
244  if (part == m_Electron) m_Tmax *= 0.5;
245  else if (part != m_Positron) {
246  G4double r = CLHEP::electron_mass_c2 / m_Mass;
247  m_Tmax = 2.0 * m_BetaGammaSq * CLHEP::electron_mass_c2 / (1.0 + 2.0 * m_Gamma * r + r * r);
248  }
249  if (m_Tmax > m_UserMaxEnergyTransfer) m_Tmax = m_UserMaxEnergyTransfer;
250  }
251  return true;
252 }
253 
254 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
255 
256 void EnergyLossForExtrapolator::Initialisation()
257 {
258  m_Initialised = true;
259  B2DEBUG(10, "EnergyLossForExtrapolator::Initialisation");
260  m_Particle = 0;
261  m_Material = 0;
262  m_KineticEnergy = 0.0;
263 
264  m_Electron = G4ParticleTable::GetParticleTable()->FindParticle("g4e_e-");
265  m_Positron = G4ParticleTable::GetParticleTable()->FindParticle("g4e_e+");
266  m_MuonPlus = G4ParticleTable::GetParticleTable()->FindParticle("g4e_mu+");
267  m_MuonMinus = G4ParticleTable::GetParticleTable()->FindParticle("g4e_mu-");
268  m_PionPlus = G4ParticleTable::GetParticleTable()->FindParticle("g4e_pi+");
269  m_PionMinus = G4ParticleTable::GetParticleTable()->FindParticle("g4e_pi-");
270  m_KaonPlus = G4ParticleTable::GetParticleTable()->FindParticle("g4e_kaon+");
271  m_KaonMinus = G4ParticleTable::GetParticleTable()->FindParticle("g4e_kaon-");
272  m_Proton = G4ParticleTable::GetParticleTable()->FindParticle("g4e_proton");
273  m_AntiProton = G4ParticleTable::GetParticleTable()->FindParticle("g4e_anti_proton");
274  m_Deuteron = G4ParticleTable::GetParticleTable()->FindParticle("g4e_deuteron");
275  m_AntiDeuteron = G4ParticleTable::GetParticleTable()->FindParticle("g4e_anti_deuteron");
276 
277  m_NMaterials = G4Material::GetNumberOfMaterials();
278  const G4MaterialTable* mtable = G4Material::GetMaterialTable();
279  m_ProductionCuts = new G4ProductionCuts();
280 
281  m_Couples.resize(m_NMaterials, 0);
282  m_Cuts.resize(m_NMaterials, DBL_MAX);
283 
284  for (G4int i = 0; i < m_NMaterials; i++) {
285  m_Couples[i] = new G4MaterialCutsCouple((*mtable)[i], m_ProductionCuts);
286  }
287 
288  m_DedxElectron = PrepareTable();
289  m_DedxPositron = PrepareTable();
290  m_DedxMuon = PrepareTable();
291  m_DedxPion = PrepareTable();
292  m_DedxKaon = PrepareTable();
293  m_DedxProton = PrepareTable();
294  m_DedxDeuteron = PrepareTable();
295  m_RangeElectron = PrepareTable();
296  m_RangePositron = PrepareTable();
297  m_RangeMuon = PrepareTable();
298  m_RangePion = PrepareTable();
299  m_RangeKaon = PrepareTable();
300  m_RangeProton = PrepareTable();
301  m_RangeDeuteron = PrepareTable();
302  m_InvRangeElectron = PrepareTable();
303  m_InvRangePositron = PrepareTable();
304  m_InvRangeMuon = PrepareTable();
305  m_InvRangePion = PrepareTable();
306  m_InvRangeKaon = PrepareTable();
307  m_InvRangeProton = PrepareTable();
308  m_InvRangeDeuteron = PrepareTable();
309  m_MscatElectron = PrepareTable();
310 
311  G4LossTableBuilder& builder = *(G4LossTableManager::Instance()->GetTableBuilder());
312 
313  B2DEBUG(10, "EnergyLossForExtrapolator Builds electron tables");
314  ComputeElectronDEDX(m_Electron, m_DedxElectron);
315  builder.BuildRangeTable(m_DedxElectron, m_RangeElectron);
316  builder.BuildInverseRangeTable(m_RangeElectron, m_InvRangeElectron);
317 
318  B2DEBUG(10, "EnergyLossForExtrapolator Builds positron tables");
319  ComputeElectronDEDX(m_Positron, m_DedxPositron);
320  builder.BuildRangeTable(m_DedxPositron, m_RangePositron);
321  builder.BuildInverseRangeTable(m_RangePositron, m_InvRangePositron);
322 
323  B2DEBUG(10, "EnergyLossForExtrapolator Builds muon tables");
324  ComputeMuonDEDX(m_MuonPlus, m_DedxMuon);
325  builder.BuildRangeTable(m_DedxMuon, m_RangeMuon);
326  builder.BuildInverseRangeTable(m_RangeMuon, m_InvRangeMuon);
327 
328  B2DEBUG(10, "EnergyLossForExtrapolator Builds pion tables");
329  ComputeHadronDEDX(m_PionPlus, m_DedxPion);
330  builder.BuildRangeTable(m_DedxPion, m_RangePion);
331  builder.BuildInverseRangeTable(m_RangePion, m_InvRangePion);
332 
333  B2DEBUG(10, "EnergyLossForExtrapolator Builds kaon tables");
334  ComputeHadronDEDX(m_KaonPlus, m_DedxKaon);
335  builder.BuildRangeTable(m_DedxKaon, m_RangeKaon);
336  builder.BuildInverseRangeTable(m_RangeKaon, m_InvRangeKaon);
337 
338  B2DEBUG(10, "EnergyLossForExtrapolator Builds proton tables");
339  ComputeHadronDEDX(m_Proton, m_DedxProton);
340  builder.BuildRangeTable(m_DedxProton, m_RangeProton);
341  builder.BuildInverseRangeTable(m_RangeProton, m_InvRangeProton);
342 
343  B2DEBUG(10, "EnergyLossForExtrapolator Builds deuteron tables");
344  ComputeHadronDEDX(m_Deuteron, m_DedxDeuteron);
345  builder.BuildRangeTable(m_DedxDeuteron, m_RangeDeuteron);
346  builder.BuildInverseRangeTable(m_RangeDeuteron, m_InvRangeDeuteron);
347 
348  ComputeTransportXS(m_Electron, m_MscatElectron);
349 }
350 
351 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
352 
353 G4PhysicsTable* EnergyLossForExtrapolator::PrepareTable()
354 {
355  G4PhysicsTable* table = new G4PhysicsTable();
356 
357  for (G4int i = 0; i < m_NMaterials; i++) {
358 
359  G4PhysicsVector* v = new G4PhysicsLogVector(m_UserTmin, m_UserTmax, m_Nbins);
360  v->SetSpline(true);
361  table->push_back(v);
362  }
363  return table;
364 }
365 
366 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
367 
368 G4double EnergyLossForExtrapolator::ComputeDEDX(G4double kinEnergy,
369  const G4ParticleDefinition* part)
370 {
371  G4double x = 0.0;
372  if (part == m_Electron) {
373  x = ComputeValue(kinEnergy, m_DedxElectron);
374  } else if (part == m_Positron) {
375  x = ComputeValue(kinEnergy, m_DedxPositron);
376  } else if (part == m_MuonPlus || part == m_MuonMinus) {
377  x = ComputeValue(kinEnergy, m_DedxMuon);
378  } else if (part == m_PionPlus || part == m_PionMinus) {
379  x = ComputeValue(kinEnergy, m_DedxPion);
380  } else if (part == m_KaonPlus || part == m_KaonMinus) {
381  x = ComputeValue(kinEnergy, m_DedxKaon);
382  } else if (part == m_Proton || part == m_AntiProton) {
383  x = ComputeValue(kinEnergy, m_DedxProton);
384  } else if (part == m_Deuteron || part == m_AntiDeuteron) {
385  x = ComputeValue(kinEnergy, m_DedxDeuteron);
386  } else {
387  B2FATAL("EnergyLossForExtrapolator::ComputeDEDX has no table for particle " << part->GetParticleName());
388  }
389  return x;
390 }
391 
392 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
393 
394 G4double EnergyLossForExtrapolator::ComputeRange(G4double kinEnergy,
395  const G4ParticleDefinition* part)
396 {
397  G4double x = 0.0;
398  if (part == m_Electron) {
399  x = ComputeValue(kinEnergy, m_RangeElectron);
400  } else if (part == m_Positron) {
401  x = ComputeValue(kinEnergy, m_RangePositron);
402  } else if (part == m_MuonPlus || part == m_MuonMinus) {
403  x = ComputeValue(kinEnergy, m_RangeMuon);
404  } else if (part == m_PionPlus || part == m_PionMinus) {
405  x = ComputeValue(kinEnergy, m_RangePion);
406  } else if (part == m_KaonPlus || part == m_KaonMinus) {
407  x = ComputeValue(kinEnergy, m_RangeKaon);
408  } else if (part == m_Proton || part == m_AntiProton) {
409  x = ComputeValue(kinEnergy, m_RangeProton);
410  } else if (part == m_Deuteron || part == m_AntiDeuteron) {
411  x = ComputeValue(kinEnergy, m_RangeDeuteron);
412  } else {
413  B2FATAL("EnergyLossForExtrapolator::ComputeRange has no table for particle " << part->GetParticleName());
414  }
415  return x;
416 }
417 
418 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
419 
420 G4double EnergyLossForExtrapolator::ComputeEnergy(G4double range,
421  const G4ParticleDefinition* part)
422 {
423  G4double x = 0.0;
424  if (part == m_Electron) {
425  x = ComputeValue(range, m_InvRangeElectron);
426  } else if (part == m_Positron) {
427  x = ComputeValue(range, m_InvRangePositron);
428  } else if (part == m_MuonPlus || part == m_MuonMinus) {
429  x = ComputeValue(range, m_InvRangeMuon);
430  } else if (part == m_PionPlus || part == m_PionMinus) {
431  x = ComputeValue(range, m_InvRangePion);
432  } else if (part == m_KaonPlus || part == m_KaonMinus) {
433  x = ComputeValue(range, m_InvRangeKaon);
434  } else if (part == m_Proton || part == m_AntiProton) {
435  x = ComputeValue(range, m_InvRangeProton);
436  } else if (part == m_Deuteron || part == m_AntiDeuteron) {
437  x = ComputeValue(range, m_InvRangeDeuteron);
438  } else {
439  B2FATAL("EnergyLossForExtrapolator::ComputeEnergy has no table for particle " << part->GetParticleName());
440  }
441  return x;
442 }
443 
444 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
445 
446 void EnergyLossForExtrapolator::ComputeElectronDEDX(const G4ParticleDefinition* part,
447  G4PhysicsTable* table)
448 {
449  G4MollerBhabhaModel* ioni = new G4MollerBhabhaModel();
450 #if G4VERSION_NUMBER < 1001
451  G4eBremsstrahlungModel* brem = new G4eBremsstrahlungModel();
452 #else
453  G4eBremsstrahlungRelModel* brem = new G4eBremsstrahlungRelModel();
454 #endif
455  G4ParticleChange* ioniPC = new G4ParticleChange();
456  ioni->SetParticleChange(ioniPC);
457  G4ParticleChange* bremPC = new G4ParticleChange();
458  brem->SetParticleChange(bremPC);
459 
460  ioni->Initialise(part, m_Cuts);
461  brem->Initialise(part, m_Cuts);
462 
463  m_Mass = CLHEP::electron_mass_c2;
464  m_ChargeSq = 1.0;
465  m_Particle = part;
466 
467  B2DEBUG(1, "EnergyLossForExtrapolator::ComputeElectronDEDX for " << part->GetParticleName());
468 
469  for (G4int i = 0; i < m_NMaterials; i++) {
470 
471  B2DEBUG(10, "EnergyLossForExtrapolator::ComputeElectronDEDX(): i= " << i << " mat= " <<
472  (*G4Material::GetMaterialTable())[i]->GetName());
473  const G4MaterialCutsCouple* couple = m_Couples[i];
474  G4PhysicsVector* aVector = (*table)[i];
475 
476  for (G4int j = 0; j <= m_Nbins; j++) {
477 
478  G4double e = aVector->Energy(j);
479  G4double dedx = ioni->ComputeDEDX(couple, part, e, e) +
480  brem->ComputeDEDX(couple, part, e, e);
481  B2DEBUG(10, "EnergyLossForExtrapolator::ComputeElectronDEDX(): j= " << j
482  << " e(MeV)= " << e / CLHEP::MeV
483  << " dedx(Mev/cm)= " << dedx * CLHEP::cm / CLHEP::MeV
484  << " dedx(Mev.cm2/g)= "
485  << dedx / ((CLHEP::MeV * (*G4Material::GetMaterialTable())[i]->GetDensity()) / (CLHEP::g / CLHEP::cm2)));
486  aVector->PutValue(j, dedx);
487  }
488  aVector->FillSecondDerivatives();
489  }
490  delete ioni; delete ioniPC;
491  delete brem; delete bremPC;
492 }
493 
494 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
495 
496 void EnergyLossForExtrapolator::ComputeMuonDEDX(const G4ParticleDefinition* part,
497  G4PhysicsTable* table)
498 {
499  G4BetheBlochModel* ioni = new G4BetheBlochModel();
500  G4MuPairProductionModel* pair = new G4MuPairProductionModel();
501  G4MuBremsstrahlungModel* brem = new G4MuBremsstrahlungModel();
502 
503  G4ParticleChange* ioniPC = new G4ParticleChange();
504  ioni->SetParticleChange(ioniPC);
505  G4ParticleChange* pairPC = new G4ParticleChange();
506  pair->SetParticleChange(pairPC);
507  G4ParticleChange* bremPC = new G4ParticleChange();
508  brem->SetParticleChange(bremPC);
509 
510  ioni->Initialise(part, m_Cuts);
511  pair->Initialise(part, m_Cuts);
512  brem->Initialise(part, m_Cuts);
513 
514  m_Mass = part->GetPDGMass();
515  m_ChargeSq = 1.0;
516  m_Particle = part;
517 
518  B2DEBUG(1, "EnergyLossForExtrapolator::ComputeMuonDEDX for " << part->GetParticleName());
519 
520  for (G4int i = 0; i < m_NMaterials; i++) {
521 
522  B2DEBUG(10, "EnergyLossForExtrapolator::ComputeMuonDEDX(): i= " << i << " mat= " <<
523  (*G4Material::GetMaterialTable())[i]->GetName());
524  const G4MaterialCutsCouple* couple = m_Couples[i];
525  G4PhysicsVector* aVector = (*table)[i];
526  for (G4int j = 0; j <= m_Nbins; j++) {
527 
528  G4double e = aVector->Energy(j);
529  G4double dedx = ioni->ComputeDEDX(couple, part, e, e) +
530  pair->ComputeDEDX(couple, part, e, e) +
531  brem->ComputeDEDX(couple, part, e, e);
532  aVector->PutValue(j, dedx);
533  B2DEBUG(10, "EnergyLossForExtrapolator::ComputeMuonDEDX(): j= " << j
534  << " e(MeV)= " << e / CLHEP::MeV
535  << " dedx(Mev/cm)= " << dedx * CLHEP::cm / CLHEP::MeV
536  << " dedx(Mev/(g/cm2)= "
537  << dedx / ((CLHEP::MeV * (*G4Material::GetMaterialTable())[i]->GetDensity()) / (CLHEP::g / CLHEP::cm2)));
538  }
539  aVector->FillSecondDerivatives();
540  }
541  delete ioni; delete ioniPC;
542  delete pair; delete pairPC;
543  delete brem; delete bremPC;
544 }
545 
546 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
547 
548 void EnergyLossForExtrapolator::ComputeHadronDEDX(const G4ParticleDefinition* part,
549  G4PhysicsTable* table)
550 {
551  G4BetheBlochModel* ioni = new G4BetheBlochModel();
552  G4ParticleChange* ioniPC = new G4ParticleChange();
553  ioni->SetParticleChange(ioniPC);
554  ioni->Initialise(part, m_Cuts);
555 
556  m_Mass = part->GetPDGMass();
557  double q = part->GetPDGCharge() / CLHEP::eplus;
558  m_ChargeSq = q * q;
559  m_Particle = part;
560 
561  B2DEBUG(1, "EnergyLossForExtrapolator::ComputeHadronDEDX for " << part->GetParticleName());
562 
563  for (G4int i = 0; i < m_NMaterials; i++) {
564 
565  B2DEBUG(10, "EnergyLossForExtrapolator::ComputeHadronDEDX(): i= " << i << " mat= " <<
566  (*G4Material::GetMaterialTable())[i]->GetName());
567  const G4MaterialCutsCouple* couple = m_Couples[i];
568  G4PhysicsVector* aVector = (*table)[i];
569  for (G4int j = 0; j <= m_Nbins; j++) {
570 
571  G4double e = aVector->Energy(j);
572  G4double dedx = ioni->ComputeDEDX(couple, part, e, e);
573  aVector->PutValue(j, dedx);
574  B2DEBUG(10, "EnergyLossForExtrapolator::ComputeHadronDEDX(): j= " << j
575  << " e(MeV)= " << e / CLHEP::MeV
576  << " dedx(Mev/cm)= " << dedx * CLHEP::cm / CLHEP::MeV
577  << " dedx(Mev.cm2/g)= "
578  << dedx / (((*G4Material::GetMaterialTable())[i]->GetDensity()) / (CLHEP::g / CLHEP::cm2)));
579  }
580  aVector->FillSecondDerivatives();
581  }
582  delete ioni; delete ioniPC;
583 }
584 
585 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
586 
587 void EnergyLossForExtrapolator::ComputeTransportXS(const G4ParticleDefinition* part,
588  G4PhysicsTable* table)
589 {
590  G4WentzelVIModel* msc = new G4WentzelVIModel();
591  G4ParticleChange* mscPC = new G4ParticleChange();
592  msc->SetParticleChange(mscPC);
593  msc->SetPolarAngleLimit(CLHEP::pi);
594  msc->Initialise(part, m_Cuts);
595 
596  m_Mass = part->GetPDGMass();
597  double q = part->GetPDGCharge() / CLHEP::eplus;
598  m_ChargeSq = q * q;
599  m_Particle = part;
600 
601  const G4MaterialTable* mtable = G4Material::GetMaterialTable();
602  B2DEBUG(1, "EnergyLossForExtrapolator::ComputeTransportXS for " << part->GetParticleName());
603 
604  for (G4int i = 0; i < m_NMaterials; i++) {
605 
606  const G4Material* mat = (*mtable)[i];
607  msc->SetCurrentCouple(m_Couples[i]);
608  B2DEBUG(10, "EnergyLossForExtrapolator::ComputeTransportXS(): i= " << i << " mat= " << mat->GetName());
609  G4PhysicsVector* aVector = (*table)[i];
610  for (G4int j = 0; j <= m_Nbins; j++) {
611 
612  G4double e = aVector->Energy(j);
613  G4double xs = msc->CrossSectionPerVolume(mat, part, e);
614  aVector->PutValue(j, xs);
615  B2DEBUG(10, "EnergyLossForExtrapolator::ComputeTransportXS(): j= " << j << " e(MeV)= " << e / CLHEP::MeV
616  << " xs(1/mm)= " << xs * CLHEP::mm);
617  }
618  aVector->FillSecondDerivatives();
619  }
620  delete msc; delete mscPC;
621 }
622 
623 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
624 
Belle2::Simulation::EnergyLossForExtrapolator::m_LinLossLimit
G4double m_LinLossLimit
Step-length limit in units of range (0.01); not modifiable by user.
Definition: EnergyLossForExtrapolator.h:295
Belle2::Simulation::EnergyLossForExtrapolator::m_KaonPlus
const G4ParticleDefinition * m_KaonPlus
Pointer to definition of the positive kaon.
Definition: EnergyLossForExtrapolator.h:169
Belle2::Simulation::EnergyLossForExtrapolator::ComputeTransportXS
void ComputeTransportXS(const G4ParticleDefinition *part, G4PhysicsTable *table)
Fill the table with the multiple-scattering cross section of a particle.
Definition: EnergyLossForExtrapolator.cc:587
Belle2::Simulation::EnergyLossForExtrapolator::ComputeMuonDEDX
void ComputeMuonDEDX(const G4ParticleDefinition *part, G4PhysicsTable *table)
Fill the table with the specific ionization energy loss of a muon.
Definition: EnergyLossForExtrapolator.cc:496
Belle2::Simulation::EnergyLossForExtrapolator::m_Particle
const G4ParticleDefinition * m_Particle
Pointer to definition of the currently cached particle.
Definition: EnergyLossForExtrapolator.h:148
Belle2::Simulation::EnergyLossForExtrapolator::ComputeDEDX
G4double ComputeDEDX(G4double kinEnergy, const G4ParticleDefinition *)
Get specific ionization energy loss for the given kinetic energy and particle.
Definition: EnergyLossForExtrapolator.cc:368
Belle2::Simulation::EnergyLossForExtrapolator::EnergyLossForExtrapolator
EnergyLossForExtrapolator()
Constructor (WITHOUT GEANT4 verbosity flag)
Definition: EnergyLossForExtrapolator.cc:47
Belle2::Simulation::EnergyLossForExtrapolator::m_DedxMuon
G4PhysicsTable * m_DedxMuon
Pointer to the muon's specific ionization energy loss vs KE table.
Definition: EnergyLossForExtrapolator.h:202
Belle2::Simulation::EnergyLossForExtrapolator::m_Mass
G4double m_Mass
Cached particle mass.
Definition: EnergyLossForExtrapolator.h:274
Belle2::Simulation::EnergyLossForExtrapolator::m_KineticEnergy
G4double m_KineticEnergy
Cached particle kinetic energy.
Definition: EnergyLossForExtrapolator.h:280
Belle2::Simulation::EnergyLossForExtrapolator::EnergyBeforeStep
G4double EnergyBeforeStep(G4double kinEnergy, G4double step, const G4Material *, const G4ParticleDefinition *)
Get kinetic energy before a step in given material by given particle.
Definition: EnergyLossForExtrapolator.cc:163
Belle2::Simulation::EnergyLossForExtrapolator::m_DedxKaon
G4PhysicsTable * m_DedxKaon
Pointer to the kaon's specific ionization energy loss vs KE table.
Definition: EnergyLossForExtrapolator.h:208
Belle2::Simulation::EnergyLossForExtrapolator::ComputeRange
G4double ComputeRange(G4double kinEnergy, const G4ParticleDefinition *)
Get range for the given kinetic energy and particle.
Definition: EnergyLossForExtrapolator.cc:394
Belle2::Simulation::EnergyLossForExtrapolator::m_DedxDeuteron
G4PhysicsTable * m_DedxDeuteron
Pointer to the deuteron's specific ionization energy loss vs KE table.
Definition: EnergyLossForExtrapolator.h:214
Belle2::Simulation::EnergyLossForExtrapolator::m_DedxProton
G4PhysicsTable * m_DedxProton
Pointer to the proton's specific ionization energy loss vs KE table.
Definition: EnergyLossForExtrapolator.h:211
Belle2::Simulation::EnergyLossForExtrapolator::m_UserTmax
G4double m_UserTmax
User's maximum kinetic energy for particles (default 10 TeV)
Definition: EnergyLossForExtrapolator.h:301
Belle2::Simulation::EnergyLossForExtrapolator::m_MscatElectron
G4PhysicsTable * m_MscatElectron
Pointer to the electron's multiple-scattering cross section vs KE table.
Definition: EnergyLossForExtrapolator.h:259
Belle2::Simulation::EnergyLossForExtrapolator::m_RangeKaon
G4PhysicsTable * m_RangeKaon
Pointer to the kaon's range vs KE table.
Definition: EnergyLossForExtrapolator.h:229
Belle2::Simulation::EnergyLossForExtrapolator::PrepareTable
G4PhysicsTable * PrepareTable()
Create a new table to store one type of kinematics data for a particle.
Definition: EnergyLossForExtrapolator.cc:353
Belle2::Simulation::EnergyLossForExtrapolator::m_Gamma
G4double m_Gamma
Cached particle's gamma value.
Definition: EnergyLossForExtrapolator.h:283
Belle2::Simulation::EnergyLossForExtrapolator::m_InvRangeProton
G4PhysicsTable * m_InvRangeProton
Pointer to the proton's inverse-range vs KE table.
Definition: EnergyLossForExtrapolator.h:253
Belle2::Simulation::EnergyLossForExtrapolator::m_ProductionCuts
G4ProductionCuts * m_ProductionCuts
Pointer to the internal cache of default G4ProductionCuts.
Definition: EnergyLossForExtrapolator.h:190
Belle2::Simulation::EnergyLossForExtrapolator::m_InvRangeKaon
G4PhysicsTable * m_InvRangeKaon
Pointer to the kaon's inverse-range vs KE table.
Definition: EnergyLossForExtrapolator.h:250
Belle2::Simulation::EnergyLossForExtrapolator::m_AntiDeuteron
const G4ParticleDefinition * m_AntiDeuteron
Pointer to definition of the antideuteron.
Definition: EnergyLossForExtrapolator.h:184
Belle2::Simulation::EnergyLossForExtrapolator::m_InvRangePositron
G4PhysicsTable * m_InvRangePositron
Pointer to the positron's inverse-range vs KE table.
Definition: EnergyLossForExtrapolator.h:241
Belle2::Simulation::EnergyLossForExtrapolator::m_RangeElectron
G4PhysicsTable * m_RangeElectron
Pointer to the electron's range vs KE table.
Definition: EnergyLossForExtrapolator.h:217
Belle2::Simulation::EnergyLossForExtrapolator::m_UserMaxEnergyTransfer
G4double m_UserMaxEnergyTransfer
User's upper limit on maximum kinetic-energy loss (default infinity)
Definition: EnergyLossForExtrapolator.h:304
Belle2::Simulation::EnergyLossForExtrapolator::m_AntiProton
const G4ParticleDefinition * m_AntiProton
Pointer to definition of the antiproton.
Definition: EnergyLossForExtrapolator.h:178
Belle2::Simulation::EnergyLossForExtrapolator::m_Deuteron
const G4ParticleDefinition * m_Deuteron
Pointer to definition of the deuteron.
Definition: EnergyLossForExtrapolator.h:181
Belle2::Simulation::EnergyLossForExtrapolator::ComputeValue
G4double ComputeValue(G4double x, const G4PhysicsTable *table)
Get the tabulated energy-loss, fluctuation, or scattering value for the given input.
Definition: EnergyLossForExtrapolator.h:403
Belle2::Simulation::EnergyLossForExtrapolator::m_KaonMinus
const G4ParticleDefinition * m_KaonMinus
Pointer to definition of the negative kaon.
Definition: EnergyLossForExtrapolator.h:172
Belle2::Simulation::EnergyLossForExtrapolator::m_MaterialIndex
G4int m_MaterialIndex
Cached material index.
Definition: EnergyLossForExtrapolator.h:265
Belle2::Simulation::EnergyLossForExtrapolator::~EnergyLossForExtrapolator
~EnergyLossForExtrapolator()
Destructor.
Definition: EnergyLossForExtrapolator.cc:108
Belle2::Simulation::EnergyLossForExtrapolator::m_DedxPositron
G4PhysicsTable * m_DedxPositron
Pointer to the positron's specific ionization energy loss vs KE table.
Definition: EnergyLossForExtrapolator.h:199
Belle2::Simulation::EnergyLossForExtrapolator::m_RadLength
G4double m_RadLength
Cached material radiation length.
Definition: EnergyLossForExtrapolator.h:271
Belle2::Simulation::EnergyLossForExtrapolator::m_BetaGammaSq
G4double m_BetaGammaSq
Cached particle's beta*gamma squared.
Definition: EnergyLossForExtrapolator.h:286
Belle2::Simulation::EnergyLossForExtrapolator::Initialisation
void Initialisation()
Initialize tables used to calculate energy loss, fluctuation, and scattering for electrons,...
Definition: EnergyLossForExtrapolator.cc:256
Belle2::Simulation::EnergyLossForExtrapolator::ComputeHadronDEDX
void ComputeHadronDEDX(const G4ParticleDefinition *part, G4PhysicsTable *table)
Fill the table with the specific ionization energy loss of a hadron.
Definition: EnergyLossForExtrapolator.cc:548
Belle2::Simulation::EnergyLossForExtrapolator::m_UserTmin
G4double m_UserTmin
User's minimum kinetic energy for particles (default 1 MeV)
Definition: EnergyLossForExtrapolator.h:298
Belle2::Simulation::EnergyLossForExtrapolator::m_InvRangeElectron
G4PhysicsTable * m_InvRangeElectron
Pointer to the electron's inverse-range vs KE table.
Definition: EnergyLossForExtrapolator.h:238
Belle2::Simulation::EnergyLossForExtrapolator::m_DedxElectron
G4PhysicsTable * m_DedxElectron
Pointer to the electron's specific ionization energy loss vs KE table.
Definition: EnergyLossForExtrapolator.h:196
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::Simulation::EnergyLossForExtrapolator::m_RangeProton
G4PhysicsTable * m_RangeProton
Pointer to the proton's range vs KE table.
Definition: EnergyLossForExtrapolator.h:232
Belle2::Simulation::EnergyLossForExtrapolator::m_MuonPlus
const G4ParticleDefinition * m_MuonPlus
Pointer to definition of the positive muon.
Definition: EnergyLossForExtrapolator.h:157
Belle2::Simulation::EnergyLossForExtrapolator::m_NMaterials
G4int m_NMaterials
Number of materials in current geometry.
Definition: EnergyLossForExtrapolator.h:310
Belle2::Simulation::EnergyLossForExtrapolator::m_PionMinus
const G4ParticleDefinition * m_PionMinus
Pointer to definition of the negative pion.
Definition: EnergyLossForExtrapolator.h:166
Belle2::Simulation::EnergyLossForExtrapolator::m_Material
const G4Material * m_Material
Pointer to internally cached material.
Definition: EnergyLossForExtrapolator.h:262
Belle2::Simulation::EnergyLossForExtrapolator::m_RangeDeuteron
G4PhysicsTable * m_RangeDeuteron
Pointer to the deuteron's range vs KE table.
Definition: EnergyLossForExtrapolator.h:235
Belle2::Simulation::EnergyLossForExtrapolator::ComputeTrueStep
G4double ComputeTrueStep(const G4Material *, const G4ParticleDefinition *part, G4double kinEnergy, G4double stepLength)
Get average scattering angle after a step in given material by given particle.
Definition: EnergyLossForExtrapolator.h:376
Belle2::Simulation::EnergyLossForExtrapolator::m_RangePion
G4PhysicsTable * m_RangePion
Pointer to the pion's range vs KE table.
Definition: EnergyLossForExtrapolator.h:226
Belle2::Simulation::EnergyLossForExtrapolator::m_PionPlus
const G4ParticleDefinition * m_PionPlus
Pointer to definition of the positive pion.
Definition: EnergyLossForExtrapolator.h:163
Belle2::Simulation::EnergyLossForExtrapolator::ComputeElectronDEDX
void ComputeElectronDEDX(const G4ParticleDefinition *part, G4PhysicsTable *table)
Fill the table with the specific ionization energy loss of an electron.
Definition: EnergyLossForExtrapolator.cc:446
Belle2::Simulation::EnergyLossForExtrapolator::m_InvRangePion
G4PhysicsTable * m_InvRangePion
Pointer to the pion's inverse-range vs KE table.
Definition: EnergyLossForExtrapolator.h:247
Belle2::Simulation::EnergyLossForExtrapolator::m_Positron
const G4ParticleDefinition * m_Positron
Pointer to definition of the positron.
Definition: EnergyLossForExtrapolator.h:154
Belle2::Simulation::EnergyLossForExtrapolator::EnergyAfterStep
G4double EnergyAfterStep(G4double kinEnergy, G4double step, const G4Material *, const G4ParticleDefinition *)
Get kinetic energy after a step in given material by given particle.
Definition: EnergyLossForExtrapolator.cc:139
Belle2::Simulation::EnergyLossForExtrapolator::ComputeEnergy
G4double ComputeEnergy(G4double range, const G4ParticleDefinition *)
Get kinetic energy corresponding to the given range and particle.
Definition: EnergyLossForExtrapolator.cc:420
Belle2::Simulation::EnergyLossForExtrapolator::m_RangeMuon
G4PhysicsTable * m_RangeMuon
Pointer to the muon's range vs KE table.
Definition: EnergyLossForExtrapolator.h:223
Belle2::Simulation::EnergyLossForExtrapolator::m_ChargeSq
G4double m_ChargeSq
Cached charge-squared (in units of e)
Definition: EnergyLossForExtrapolator.h:277
Belle2::Simulation::EnergyLossForExtrapolator::m_Couples
std::vector< const G4MaterialCutsCouple * > m_Couples
List of material-cuts pairings.
Definition: EnergyLossForExtrapolator.h:193
Belle2::Simulation::EnergyLossForExtrapolator::m_InvRangeMuon
G4PhysicsTable * m_InvRangeMuon
Pointer to the muon's inverse-range vs KE table.
Definition: EnergyLossForExtrapolator.h:244
Belle2::Simulation::EnergyLossForExtrapolator::m_InvRangeDeuteron
G4PhysicsTable * m_InvRangeDeuteron
Pointer to the deuteron's inverse-range vs KE table.
Definition: EnergyLossForExtrapolator.h:256
Belle2::Simulation::EnergyLossForExtrapolator::m_RangePositron
G4PhysicsTable * m_RangePositron
Pointer to the positron's range vs KE table.
Definition: EnergyLossForExtrapolator.h:220
Belle2::Simulation::EnergyLossForExtrapolator::m_ElectronDensity
G4double m_ElectronDensity
Cached material electron density.
Definition: EnergyLossForExtrapolator.h:268
Belle2::Simulation::EnergyLossForExtrapolator::m_BetaSq
G4double m_BetaSq
Cached particle's beta squared.
Definition: EnergyLossForExtrapolator.h:289
Belle2::Simulation::EnergyLossForExtrapolator::m_Electron
const G4ParticleDefinition * m_Electron
Pointer to definition of the electron.
Definition: EnergyLossForExtrapolator.h:151
Belle2::Simulation::EnergyLossForExtrapolator::SetupKinematics
G4bool SetupKinematics(const G4ParticleDefinition *, const G4Material *, G4double kinEnergy)
Save current particle properties, kinetic energy and material in internal cached state.
Definition: EnergyLossForExtrapolator.cc:210
Belle2::Simulation::EnergyLossForExtrapolator::TrueStepLength
G4double TrueStepLength(G4double kinEnergy, G4double step, const G4Material *, const G4ParticleDefinition *part)
Get true step length for given particle and kinetic energy.
Definition: EnergyLossForExtrapolator.cc:187
Belle2::Simulation::EnergyLossForExtrapolator::m_Proton
const G4ParticleDefinition * m_Proton
Pointer to definition of the proton.
Definition: EnergyLossForExtrapolator.h:175
Belle2::Simulation::EnergyLossForExtrapolator::m_Nbins
G4int m_Nbins
Number of bins in each energy loss, fluctuation, and multiple-scattering table (70,...
Definition: EnergyLossForExtrapolator.h:307
Belle2::Simulation::EnergyLossForExtrapolator::m_Cuts
G4DataVector m_Cuts
Vector of particle cuts.
Definition: EnergyLossForExtrapolator.h:187
Belle2::Simulation::EnergyLossForExtrapolator::m_MuonMinus
const G4ParticleDefinition * m_MuonMinus
Pointer to definition of the negative muon.
Definition: EnergyLossForExtrapolator.h:160
Belle2::Simulation::EnergyLossForExtrapolator::m_Tmax
G4double m_Tmax
Cached particle's maximum kinetic-energy loss.
Definition: EnergyLossForExtrapolator.h:292
Belle2::Simulation::EnergyLossForExtrapolator::m_Initialised
G4bool m_Initialised
UNUSED Flag to indicate that Initialisation() method has been called.
Definition: EnergyLossForExtrapolator.h:313
Belle2::Simulation::EnergyLossForExtrapolator::m_DedxPion
G4PhysicsTable * m_DedxPion
Pointer to the pion's specific ionization energy loss vs KE table.
Definition: EnergyLossForExtrapolator.h:205