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