LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
NestAlg.cxx
Go to the documentation of this file.
1 #define AVO 6.022e23 //Avogadro's number (#/mol)
2 #define EMASS 9.109e-31 * CLHEP::kg
3 #define MillerDriftSpeed true
4 
5 #define GASGAP 0.25 * CLHEP::cm //S2 generation region
6 #define BORDER 0 * CLHEP::cm //liquid-gas border z-coordinate
7 
8 #define QE_EFF 1 //a base or maximum quantum efficiency
9 #define phe_per_e 1 //S2 gain for quick studies
10 
11 // different field regions, for gamma-X studies
12 #define WIN 0 * CLHEP::mm //top Cu block (also, quartz window)
13 #define TOP 0 //top grid wires
14 #define ANE 0 //anode mesh
15 #define SRF 0 //liquid-gas interface
16 #define GAT 0 //gate grid
17 #define CTH 0 //cathode grid
18 #define BOT 0 //bottom PMT grid
19 #define PMT 0 //bottom Cu block and PMTs
20 #define MIN_ENE -1 * CLHEP::eV //lets you turn NEST off BELOW a certain energy
21 #define MAX_ENE 1. * CLHEP::TeV //lets you turn NEST off ABOVE a certain energy
22 #define HIENLIM 5 * CLHEP::MeV //energy at which Doke model used exclusively
23 #define R_TOL 0.2 * CLHEP::mm //tolerance (for edge events)
24 #define R_MAX 1 * CLHEP::km //for corraling diffusing electrons
25 #define Density_LXe 2.9 //reference density for density-dep. effects
26 #define Density_LAr 1.393
27 #define Density_LNe 1.207
28 #define Density_LKr 2.413
29 
30 #include "Geant4/G4Ions.hh"
31 #include "Geant4/G4OpticalPhoton.hh"
32 #include "Geant4/G4VProcess.hh"
33 
36 
37 #include "CLHEP/Random/RandFlat.h"
38 #include "CLHEP/Random/RandGauss.h"
39 
40 G4bool diffusion = true;
41 
42 G4bool SinglePhase = false, ThomasImelTail = true, OutElectrons = true;
43 
44 G4double biExc = 0.77; //for alpha particles (bi-excitonic collisions)
45 
46 //----------------------------------------------------------------------------
47 // Default constructor will return no photons or electrons unless the set
48 // methods are called.
49 NestAlg::NestAlg(CLHEP::HepRandomEngine& engine)
50  : fYieldFactor(0)
51  , fExcitationRatio(0)
52  , fNumScintPhotons(0)
53  , fNumIonElectrons(0)
54  , fEnergyDep(0.)
55  , fEngine(engine)
56 {
57  fElementPropInit[2] = false;
58  fElementPropInit[10] = false;
59  fElementPropInit[18] = false;
60  fElementPropInit[36] = false;
61  fElementPropInit[54] = false;
62 }
63 
64 //----------------------------------------------------------------------------
65 NestAlg::NestAlg(double yieldFactor, CLHEP::HepRandomEngine& engine)
66  : fYieldFactor(yieldFactor)
67  , fExcitationRatio(0.)
68  , fNumScintPhotons(0)
69  , fNumIonElectrons(0)
70  , fEnergyDep(0.)
71  , fEngine(engine)
72 {
73  fElementPropInit[2] = false;
74  fElementPropInit[10] = false;
75  fElementPropInit[18] = false;
76  fElementPropInit[36] = false;
77  fElementPropInit[54] = false;
78 }
79 
80 //----------------------------------------------------------------------------
81 const G4VParticleChange& NestAlg::CalculateIonizationAndScintillation(G4Track const& aTrack,
82  G4Step const& aStep)
83 {
84  CLHEP::RandGauss GaussGen(fEngine);
85  CLHEP::RandFlat UniformGen(fEngine);
86 
87  // reset the variables accessed by other objects
88  // make the energy deposit the energy in this step,
89  // set the number of electrons and photons to 0
90  fEnergyDep = aStep.GetTotalEnergyDeposit();
91  fNumIonElectrons = 0.;
92  fNumScintPhotons = 0.;
93 
94  if (!fYieldFactor) //set YF=0 when you want S1Light off in your sim
95  return fParticleChange;
96 
97  bool fTrackSecondariesFirst = false;
98  bool fExcitedNucleus = false;
99  bool fVeryHighEnergy = false;
100  bool fAlpha = false;
101  bool fMultipleScattering = false;
102  double fKr83m = 0.;
103 
104  if (aTrack.GetParentID() == 0 && aTrack.GetCurrentStepNumber() == 1) {
105  fExcitedNucleus = false; //an initialization or reset
106  fVeryHighEnergy = false; //initializes or (later) resets this
107  fAlpha = false; //ditto
108  fMultipleScattering = false;
109  }
110 
111  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
112  G4ParticleDefinition* pDef = aParticle->GetDefinition();
113  G4String particleName = pDef->GetParticleName();
114  const G4Material* aMaterial = aStep.GetPreStepPoint()->GetMaterial();
115  const G4Material* bMaterial = aStep.GetPostStepPoint()->GetMaterial();
116 
117  if ((particleName == "neutron" || particleName == "antineutron") &&
118  aStep.GetTotalEnergyDeposit() <= 0)
119  return fParticleChange;
120 
121  // code for determining whether the present/next material is noble
122  // element, or, in other words, for checking if either is a valid NEST
123  // scintillating material, and save Z for later L calculation, or
124  // return if no valid scintillators are found on this step, which is
125  // protection against G4Exception or seg. fault/violation
126  G4Element *ElementA = NULL, *ElementB = NULL;
127  if (aMaterial) {
128  const G4ElementVector* theElementVector1 = aMaterial->GetElementVector();
129  ElementA = (*theElementVector1)[0];
130  }
131  if (bMaterial) {
132  const G4ElementVector* theElementVector2 = bMaterial->GetElementVector();
133  ElementB = (*theElementVector2)[0];
134  }
135  G4int z1, z2, j = 1;
136  G4bool NobleNow = false, NobleLater = false;
137  if (ElementA)
138  z1 = (G4int)(ElementA->GetZ());
139  else
140  z1 = -1;
141  if (ElementB)
142  z2 = (G4int)(ElementB->GetZ());
143  else
144  z2 = -1;
145  if (z1 == 2 || z1 == 10 || z1 == 18 || z1 == 36 || z1 == 54) {
146  NobleNow = true;
147  // j = (G4int)aMaterial->GetMaterialPropertiesTable()->
148  // GetConstProperty("TOTALNUM_INT_SITES"); //get current number
149  //if ( j < 0 ) {
150  if (aTrack.GetParentID() == 0 && !fElementPropInit[z1]) {
151  InitMatPropValues(aMaterial->GetMaterialPropertiesTable(), z1);
152  j = 0; //no sites yet
153  } //material properties initialized
154  } //end of atomic number check
155  if (z2 == 2 || z2 == 10 || z2 == 18 || z2 == 36 || z2 == 54) {
156  NobleLater = true;
157  // j = (G4int)bMaterial->GetMaterialPropertiesTable()->
158  // GetConstProperty("TOTALNUM_INT_SITES");
159  //if ( j < 0 ) {
160  if (aTrack.GetParentID() == 0 && !fElementPropInit[z2]) {
161  InitMatPropValues(bMaterial->GetMaterialPropertiesTable(), z2);
162  j = 0; //no sites yet
163  } //material properties initialized
164  } //end of atomic number check
165 
166  if (!NobleNow && !NobleLater) return fParticleChange;
167 
168  // retrieval of the particle's position, time, attributes at both the
169  // beginning and the end of the current step along its track
170  G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
171  G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
172  G4ThreeVector x1 = pPostStepPoint->GetPosition();
173  G4ThreeVector x0 = pPreStepPoint->GetPosition();
174  G4double evtStrt = pPreStepPoint->GetGlobalTime();
175  G4double t0 = pPreStepPoint->GetLocalTime();
176  G4double t1 = pPostStepPoint->GetLocalTime();
177 
178  // now check if we're entering a scintillating material (inside) or
179  // leaving one (outside), in order to determine (later on in the code,
180  // based on the booleans inside & outside) whether to add/subtract
181  // energy that can potentially be deposited from the system
182  G4bool outside = false, inside = false, InsAndOuts = false;
183  G4MaterialPropertiesTable* aMaterialPropertiesTable = aMaterial->GetMaterialPropertiesTable();
184  if (NobleNow && !NobleLater) outside = true;
185  if (!NobleNow && NobleLater) {
186  aMaterial = bMaterial;
187  inside = true;
188  z1 = z2;
189  aMaterialPropertiesTable = bMaterial->GetMaterialPropertiesTable();
190  }
191  if (NobleNow && NobleLater && aMaterial->GetDensity() != bMaterial->GetDensity())
192  InsAndOuts = true;
193 
194  // retrieve scintillation-related material properties
195  G4double Density = aMaterial->GetDensity() / (CLHEP::g / CLHEP::cm3);
196  G4double nDensity = Density * AVO; //molar mass factor applied below
197  G4int Phase = aMaterial->GetState(); //solid, liquid, or gas?
198  G4double ElectricField(0.), FieldSign(0.); //for field quenching of S1
199  G4bool GlobalFields = false;
200 
201  if ((WIN == 0) && !TOP && !ANE && !SRF && !GAT && !CTH && !BOT && !PMT) {
202  ElectricField = aMaterialPropertiesTable->GetConstProperty("ELECTRICFIELD");
203  GlobalFields = true;
204  }
205  else {
206  if (x1[2] < WIN && x1[2] > TOP && Phase == kStateGas)
207  ElectricField = aMaterialPropertiesTable->GetConstProperty("ELECTRICFIELDWINDOW");
208  else if (x1[2] < TOP && x1[2] > ANE && Phase == kStateGas)
209  ElectricField = aMaterialPropertiesTable->GetConstProperty("ELECTRICFIELDTOP");
210  else if (x1[2] < ANE && x1[2] > SRF && Phase == kStateGas)
211  ElectricField = aMaterialPropertiesTable->GetConstProperty("ELECTRICFIELDANODE");
212  else if (x1[2] < SRF && x1[2] > GAT && Phase == kStateLiquid)
213  ElectricField = aMaterialPropertiesTable->GetConstProperty("ELECTRICFIELDSURFACE");
214  else if (x1[2] < GAT && x1[2] > CTH && Phase == kStateLiquid)
215  ElectricField = aMaterialPropertiesTable->GetConstProperty("ELECTRICFIELDGATE");
216  else if (x1[2] < CTH && x1[2] > BOT && Phase == kStateLiquid)
217  ElectricField = aMaterialPropertiesTable->GetConstProperty("ELECTRICFIELDCATHODE");
218  else if (x1[2] < BOT && x1[2] > PMT && Phase == kStateLiquid)
219  ElectricField = aMaterialPropertiesTable->GetConstProperty("ELECTRICFIELDBOTTOM");
220  else
221  ElectricField = aMaterialPropertiesTable->GetConstProperty("ELECTRICFIELD");
222  }
223  if (ElectricField >= 0)
224  FieldSign = 1;
225  else
226  FieldSign = -1;
227  ElectricField = std::abs((1e3 * ElectricField) / (CLHEP::kilovolt / CLHEP::cm));
228  G4double Temperature = aMaterial->GetTemperature();
229  G4double ScintillationYield, ResolutionScale, R0 = 1.0 * CLHEP::um, DokeBirks[3],
230  ThomasImel = 0.00, delta = 1 * CLHEP::mm;
231  DokeBirks[0] = 0.00;
232  DokeBirks[2] = 1.00;
233  G4double PhotMean = 7 * CLHEP::eV, PhotWidth = 1.0 * CLHEP::eV; //photon properties
234  G4double SingTripRatioR, SingTripRatioX, tau1, tau3, tauR = 0 * CLHEP::ns;
235  switch (z1) { //sort prop. by noble element atomic#
236  case 2: //helium
237  ScintillationYield = 1 / (41.3 * CLHEP::eV); //all W's from noble gas book
238  fExcitationRatio = 0.00; //nominal (true value unknown)
239  ResolutionScale = 0.2; //Aprile, Bolotnikov, Bolozdynya, Doke
240  PhotMean = 15.9 * CLHEP::eV;
241  tau1 = GaussGen.fire(10.0 * CLHEP::ns, 0.0 * CLHEP::ns);
242  tau3 = 1.6e3 * CLHEP::ns;
243  tauR = GaussGen.fire(13.00 * CLHEP::s, 2.00 * CLHEP::s); //McKinsey et al. 2003
244  break;
245  case 10: //neon
246  ScintillationYield = 1 / (29.2 * CLHEP::eV);
247  fExcitationRatio = 0.00; //nominal (true value unknown)
248  ResolutionScale = 0.13; //Aprile et. al book
249  PhotMean = 15.5 * CLHEP::eV;
250  PhotWidth = 0.26 * CLHEP::eV;
251  tau1 = GaussGen.fire(10.0 * CLHEP::ns, 10. * CLHEP::ns);
252  tau3 = GaussGen.fire(15.4e3 * CLHEP::ns, 200 * CLHEP::ns); //Nikkel et al. 2008
253  break;
254  case 18: //argon
255  ScintillationYield = 1 / (19.5 * CLHEP::eV);
256  fExcitationRatio = 0.21; //Aprile et. al book
257  ResolutionScale = 0.107; //Doke 1976
258  R0 = 1.568 * CLHEP::um; //Mozumder 1995
259  if (ElectricField) {
260  ThomasImel = 0.156977 * pow(ElectricField, -0.1);
261  DokeBirks[0] = 0.07 * pow((ElectricField / 1.0e3), -0.85);
262  DokeBirks[2] = 0.00;
263  }
264  else {
265  ThomasImel = 0.099;
266  DokeBirks[0] = 0.0003;
267  DokeBirks[2] = 0.75;
268  }
269  nDensity /= 39.948; //molar mass in grams per mole
270  PhotMean = 9.69 * CLHEP::eV;
271  PhotWidth = 0.22 * CLHEP::eV;
272  tau1 = GaussGen.fire(6.5 * CLHEP::ns, 0.8 * CLHEP::ns); //err from wgted avg.
273  tau3 = GaussGen.fire(1300 * CLHEP::ns, 50 * CLHEP::ns); //ibid.
274  tauR = GaussGen.fire(0.8 * CLHEP::ns, 0.2 * CLHEP::ns); //Kubota 1979
275  biExc = 0.6;
276  break;
277  case 36: //krypton
278  if (Phase == kStateGas)
279  ScintillationYield = 1 / (30.0 * CLHEP::eV);
280  else
281  ScintillationYield = 1 / (15.0 * CLHEP::eV);
282  fExcitationRatio = 0.08; //Aprile et. al book
283  ResolutionScale = 0.05; //Doke 1976
284  PhotMean = 8.43 * CLHEP::eV;
285  tau1 = GaussGen.fire(2.8 * CLHEP::ns, .04 * CLHEP::ns);
286  tau3 = GaussGen.fire(93. * CLHEP::ns, 1.1 * CLHEP::ns);
287  tauR = GaussGen.fire(12. * CLHEP::ns, .76 * CLHEP::ns);
288  break;
289  case 54: //xenon
290  default:
291  nDensity /= 131.293;
292  ScintillationYield = 48.814 + 0.80877 * Density + 2.6721 * pow(Density, 2.);
293  ScintillationYield /= CLHEP::keV; //Units: [#quanta(ph/e-) per keV]
294  //W = 13.7 eV for all recoils, Dahl thesis (that's @2.84 g/cm^3)
295  //the exciton-to-ion ratio (~0.06 for LXe at 3 g/cm^3)
296  fExcitationRatio = 0.4 - 0.11131 * Density - 0.0026651 * pow(Density, 2.);
297  ResolutionScale = 1.00 * //Fano factor <<1
298  (0.12724 - 0.032152 * Density - 0.0013492 * pow(Density, 2.));
299  //~0.1 for GXe w/ formula from Bolotnikov et al. 1995
300  if (Phase == kStateLiquid) {
301  ResolutionScale *= 1.5; //to get it to be ~0.03 for LXe
302  R0 = 16.6 * CLHEP::um; //for zero electric field
303  //length scale above which Doke model used instead of Thomas-Imel
304  if (ElectricField) //change it with field (see NEST paper)
305  R0 = 69.492 * pow(ElectricField, -0.50422) * CLHEP::um;
306  if (ElectricField) { //formulae & values all from NEST paper
307  DokeBirks[0] = 19.171 * pow(ElectricField + 25.552, -0.83057) + 0.026772;
308  DokeBirks[2] = 0.00; //only volume recombination (above)
309  }
310  else { //zero electric field magnitude
311  DokeBirks[0] = 0.18; //volume/columnar recombination factor (A)
312  DokeBirks[2] = 0.58; //geminate/Onsager recombination (C)
313  }
314  //"ThomasImel" is alpha/(a^2*v), the recomb. coeff.
315  ThomasImel = 0.05; //aka xi/(4*N_i) from the NEST paper
316  //distance used to determine when one is at a new interaction site
317  delta = 0.4 * CLHEP::mm; //distance ~30 keV x-ray travels in LXe
318  PhotMean = 6.97 * CLHEP::eV;
319  PhotWidth = 0.23 * CLHEP::eV;
320  // 178+/-14nmFWHM, taken from Jortner JchPh 42 '65.
321  //these singlet and triplet times may not be the ones you're
322  //used to, but are the world average: Kubota 79, Hitachi 83 (2
323  //data sets), Teymourian 11, Morikawa 89, and Akimov '02
324  tau1 = GaussGen.fire(3.1 * CLHEP::ns, .7 * CLHEP::ns); //err from wgted avg.
325  tau3 = GaussGen.fire(24. * CLHEP::ns, 1. * CLHEP::ns); //ibid.
326  } //end liquid
327  else if (Phase == kStateGas) {
328  if (!fAlpha)
329  fExcitationRatio = 0.07; //Nygren NIM A 603 (2009) p. 340
330  else {
331  biExc = 1.00;
332  ScintillationYield = 1 / (12.98 * CLHEP::eV);
333  } //Saito 2003
334  R0 = 0.0 * CLHEP::um; //all Doke/Birks interactions (except for alphas)
335  G4double Townsend = (ElectricField / nDensity) * 1e17;
336  DokeBirks[0] = 0.0000; //all geminate (except at zero, low fields)
337  DokeBirks[2] = 0.1933 * pow(Density, 2.6199) + 0.29754 -
338  (0.045439 * pow(Density, 2.4689) + 0.066034) * log10(ElectricField);
339  if (ElectricField > 6990) DokeBirks[2] = 0.0;
340  if (ElectricField < 1000) DokeBirks[2] = 0.2;
341  if (ElectricField < 100.) {
342  DokeBirks[0] = 0.18;
343  DokeBirks[2] = 0.58;
344  }
345  if (Density < 0.061)
346  ThomasImel = 0.041973 * pow(Density, 1.8105);
347  else if (Density >= 0.061 && Density <= 0.167)
348  ThomasImel = 5.9583e-5 + 0.0048523 * Density - 0.023109 * pow(Density, 2.);
349  else
350  ThomasImel = 6.2552e-6 * pow(Density, -1.9963);
351  if (ElectricField) ThomasImel = 1.2733e-5 * pow(Townsend / Density, -0.68426);
352  // field\density dependence from Kobayashi 2004 and Saito 2003
353  PhotMean = 7.1 * CLHEP::eV;
354  PhotWidth = 0.2 * CLHEP::eV;
355  tau1 = GaussGen.fire(5.18 * CLHEP::ns, 1.55 * CLHEP::ns);
356  tau3 = GaussGen.fire(100.1 * CLHEP::ns, 7.9 * CLHEP::ns);
357  } //end gas information (preliminary guesses)
358  else {
359  tau1 = 3.5 * CLHEP::ns;
360  tau3 = 20. * CLHEP::ns;
361  tauR = 40. * CLHEP::ns;
362  } //solid Xe
363  }
364 
365  // log present and running tally of energy deposition in this section
366  G4double anExcitationEnergy =
367  ((const G4Ions*)(pDef))->GetExcitationEnergy(); //grab nuclear energy level
368  G4double TotalEnergyDeposit = //total energy deposited so far
369  aMaterialPropertiesTable->GetConstProperty("ENERGY_DEPOSIT_TOT");
370  G4bool convert = false, annihil = false;
371  //set up special cases for pair production and positron annihilation
372  if (pPreStepPoint->GetKineticEnergy() >= (2 * CLHEP::electron_mass_c2) &&
373  !pPostStepPoint->GetKineticEnergy() && !aStep.GetTotalEnergyDeposit() &&
374  aParticle->GetPDGcode() == 22) {
375  convert = true;
376  TotalEnergyDeposit = CLHEP::electron_mass_c2;
377  }
378  if (pPreStepPoint->GetKineticEnergy() && !pPostStepPoint->GetKineticEnergy() &&
379  aParticle->GetPDGcode() == -11) {
380  annihil = true;
381  TotalEnergyDeposit += aStep.GetTotalEnergyDeposit();
382  }
383  G4bool either = false;
384  if (inside || outside || convert || annihil || InsAndOuts) either = true;
385  //conditions for returning when energy deposits too low
386  if (anExcitationEnergy < 100 * CLHEP::eV && aStep.GetTotalEnergyDeposit() < 1 * CLHEP::eV &&
387  !either && !fExcitedNucleus)
388  return fParticleChange;
389  //add current deposit to total energy budget
390  if (!annihil) TotalEnergyDeposit += aStep.GetTotalEnergyDeposit();
391  if (!convert)
392  aMaterialPropertiesTable->AddConstProperty("ENERGY_DEPOSIT_TOT", TotalEnergyDeposit);
393  //save current deposit for determining number of quanta produced now
394  TotalEnergyDeposit = aStep.GetTotalEnergyDeposit();
395 
396  // check what the current "goal" E is for dumping scintillation,
397  // often the initial kinetic energy of the parent particle, and deal
398  // with all other energy-related matters in this block of code
399  G4double InitialKinetEnergy = aMaterialPropertiesTable->GetConstProperty("ENERGY_DEPOSIT_GOL");
400  //if zero, add up initial potential and kinetic energies now
401  if (InitialKinetEnergy == 0) {
402  G4double tE = pPreStepPoint->GetKineticEnergy() + anExcitationEnergy;
403  if ((fabs(tE - 1.8 * CLHEP::keV) < CLHEP::eV || fabs(tE - 9.4 * CLHEP::keV) < CLHEP::eV) &&
404  Phase == kStateLiquid && z1 == 54)
405  tE = 9.4 * CLHEP::keV;
406  if (fKr83m && ElectricField != 0) DokeBirks[2] = 0.10;
407  aMaterialPropertiesTable->AddConstProperty("ENERGY_DEPOSIT_GOL", tE);
408  //excited nucleus is special case where accuracy reduced for total
409  //energy deposition because of G4 inaccuracies and scintillation is
410  //forced-dumped when that nucleus is fully de-excited
411  if (anExcitationEnergy) fExcitedNucleus = true;
412  }
413  //if a particle is leaving, remove its kinetic energy from the goal
414  //energy, as this will never get deposited (if depositable)
415  if (outside) {
416  aMaterialPropertiesTable->AddConstProperty(
417  "ENERGY_DEPOSIT_GOL", InitialKinetEnergy - pPostStepPoint->GetKineticEnergy());
418  if (aMaterialPropertiesTable->GetConstProperty("ENERGY_DEPOSIT_GOL") < 0)
419  aMaterialPropertiesTable->AddConstProperty("ENERGY_DEPOSIT_GOL", 0);
420  }
421  //if a particle is coming back into your scintillator, then add its
422  //energy to the goal energy
423  if (inside) {
424  aMaterialPropertiesTable->AddConstProperty(
425  "ENERGY_DEPOSIT_GOL", InitialKinetEnergy + pPreStepPoint->GetKineticEnergy());
426  if (TotalEnergyDeposit > 0 && InitialKinetEnergy == 0) {
427  aMaterialPropertiesTable->AddConstProperty("ENERGY_DEPOSIT_GOL", 0);
428  TotalEnergyDeposit = .000000;
429  }
430  }
431  if (InsAndOuts) {
432  //G4double dribble = pPostStepPoint->GetKineticEnergy() -
433  //pPreStepPoint->GetKineticEnergy();
434  aMaterialPropertiesTable->AddConstProperty("ENERGY_DEPOSIT_GOL",
435  (-0.1 * CLHEP::keV) + InitialKinetEnergy -
436  pPostStepPoint->GetKineticEnergy());
437  InitialKinetEnergy =
438  bMaterial->GetMaterialPropertiesTable()->GetConstProperty("ENERGY_DEPOSIT_GOL");
439  bMaterial->GetMaterialPropertiesTable()->AddConstProperty(
440  "ENERGY_DEPOSIT_GOL",
441  (-0.1 * CLHEP::keV) + InitialKinetEnergy + pPreStepPoint->GetKineticEnergy());
442  if (aMaterialPropertiesTable->GetConstProperty("ENERGY_DEPOSIT_GOL") < 0)
443  aMaterialPropertiesTable->AddConstProperty("ENERGY_DEPOSIT_GOL", 0);
444  if (bMaterial->GetMaterialPropertiesTable()->GetConstProperty("ENERGY_DEPOSIT_GOL") < 0)
445  bMaterial->GetMaterialPropertiesTable()->AddConstProperty("ENERGY_DEPOSIT_GOL", 0);
446  }
447  InitialKinetEnergy =
448  aMaterialPropertiesTable->GetConstProperty("ENERGY_DEPOSIT_GOL"); //grab current goal E
449  if (annihil) //if an annihilation occurred, add energy of two gammas
450  InitialKinetEnergy += 2 * CLHEP::electron_mass_c2;
451  //if pair production occurs, then subtract energy to cancel with the
452  //energy that will be added in the line above when the e+ dies
453  if (convert) InitialKinetEnergy -= 2 * CLHEP::electron_mass_c2;
454  //update the relevant material property (goal energy)
455  aMaterialPropertiesTable->AddConstProperty("ENERGY_DEPOSIT_GOL", InitialKinetEnergy);
456  if (anExcitationEnergy < 1e-100 && aStep.GetTotalEnergyDeposit() == 0 &&
457  aMaterialPropertiesTable->GetConstProperty("ENERGY_DEPOSIT_GOL") == 0 &&
458  aMaterialPropertiesTable->GetConstProperty("ENERGY_DEPOSIT_TOT") == 0)
459  return fParticleChange;
460 
461  G4String procName;
462  if (aTrack.GetCreatorProcess())
463  procName = aTrack.GetCreatorProcess()->GetProcessName();
464  else
465  procName = "NULL";
466  if (procName == "eBrem" && outside && !OutElectrons) fMultipleScattering = true;
467 
468  // next 2 codeblocks deal with position-related things
469  if (fAlpha) delta = 1000. * CLHEP::km;
470  G4int i, k, counter = 0;
471  G4double pos[3];
472  if (outside) { //leaving
473  if (aParticle->GetPDGcode() == 11 && !OutElectrons) fMultipleScattering = true;
474  x1 = x0; //prevents generation of quanta outside active volume
475  } //no scint. for e-'s that leave
476 
477  char xCoord[80];
478  char yCoord[80];
479  char zCoord[80];
480  G4bool exists = false; //for querying whether set-up of new site needed
481  for (i = 0; i < j; i++) { //loop over all saved interaction sites
482  counter = i; //save site# for later use in storing properties
483  sprintf(xCoord, "POS_X_%d", i);
484  sprintf(yCoord, "POS_Y_%d", i);
485  sprintf(zCoord, "POS_Z_%d", i);
486  pos[0] = aMaterialPropertiesTable->GetConstProperty(xCoord);
487  pos[1] = aMaterialPropertiesTable->GetConstProperty(yCoord);
488  pos[2] = aMaterialPropertiesTable->GetConstProperty(zCoord);
489  if (sqrt(pow(x1[0] - pos[0], 2.) + pow(x1[1] - pos[1], 2.) + pow(x1[2] - pos[2], 2.)) < delta) {
490  exists = true;
491  break; //we find interaction is close to an old one
492  }
493  }
494  if (!exists && TotalEnergyDeposit) { //current interaction too far away
495  counter = j;
496  sprintf(xCoord, "POS_X_%i", j);
497  sprintf(yCoord, "POS_Y_%i", j);
498  sprintf(zCoord, "POS_Z_%i", j);
499  //save 3-space coordinates of the new interaction site
500  aMaterialPropertiesTable->AddConstProperty(xCoord, x1[0]);
501  aMaterialPropertiesTable->AddConstProperty(yCoord, x1[1]);
502  aMaterialPropertiesTable->AddConstProperty(zCoord, x1[2]);
503  j++; //increment number of sites
504  aMaterialPropertiesTable-> //save
505  AddConstProperty("TOTALNUM_INT_SITES", j);
506  }
507 
508  // this is where nuclear recoil "L" factor is handled: total yield is
509  // reduced for nuclear recoil as per Lindhard theory
510 
511  //we assume you have a mono-elemental scintillator only
512  //now, grab A's and Z's of current particle and of material (avg)
513  G4double a1 = ElementA->GetA();
514  z2 = pDef->GetAtomicNumber();
515  G4double a2 = (G4double)(pDef->GetAtomicMass());
516  if (particleName == "alpha" || (z2 == 2 && a2 == 4))
517  fAlpha = true; //used later to get S1 pulse shape correct for alpha
518  if (fAlpha || abs(aParticle->GetPDGcode()) == 2112) a2 = a1; //get average A for element at hand
519  G4double epsilon = 11.5 * (TotalEnergyDeposit / CLHEP::keV) * pow(z1, (-7. / 3.));
520  G4double gamma = 3. * pow(epsilon, 0.15) + 0.7 * pow(epsilon, 0.6) + epsilon;
521  G4double kappa = 0.133 * pow(z1, (2. / 3.)) * pow(a2, (-1. / 2.)) * (2. / 3.);
522  //check if we are dealing with nuclear recoil (Z same as material)
523  if ((z1 == z2 && pDef->GetParticleType() == "nucleus") || particleName == "neutron" ||
524  particleName == "antineutron") {
525  fYieldFactor = (kappa * gamma) / (1 + kappa * gamma); //Lindhard factor
526  if (z1 == 18 && Phase == kStateLiquid)
527  fYieldFactor = 0.23 * (1 + exp(-5 * epsilon)); //liquid argon L_eff
528  //just a few safety checks, like for recombProb below
529  if (fYieldFactor > 1) fYieldFactor = 1;
530  if (fYieldFactor < 0) fYieldFactor = 0;
531  if (ElectricField == 0 && Phase == kStateLiquid) {
532  if (z1 == 54) ThomasImel = 0.19;
533  if (z1 == 18) ThomasImel = 0.25;
534  } //special TIB parameters for nuclear recoil only, in LXe and LAr
535  fExcitationRatio = 0.69337 + 0.3065 * exp(-0.008806 * pow(ElectricField, 0.76313));
536  }
537  else
538  fYieldFactor = 1.000; //default
539 
540  // determine ultimate #quanta from current E-deposition (ph+e-)
541  G4double MeanNumberOfQuanta = //total mean number of exc/ions
542  ScintillationYield * TotalEnergyDeposit;
543  //the total number of either quanta produced is equal to product of the
544  //work function, the energy deposited, and yield reduction, for NR
545  G4double sigma = sqrt(ResolutionScale * MeanNumberOfQuanta); //Fano
546  G4int NumQuanta = //stochastic variation in NumQuanta
547  G4int(floor(GaussGen.fire(MeanNumberOfQuanta, sigma) + 0.5));
548  G4double LeffVar = GaussGen.fire(fYieldFactor, 0.25 * fYieldFactor);
549  if (LeffVar > 1) { LeffVar = 1.00000; }
550  if (LeffVar < 0) { LeffVar = 0; }
551  if (fYieldFactor < 1) NumQuanta = BinomFluct(NumQuanta, LeffVar);
552  //if E below work function, can't make any quanta, and if NumQuanta
553  //less than zero because Gaussian fluctuated low, update to zero
554  if (TotalEnergyDeposit < 1 / ScintillationYield || NumQuanta < 0) NumQuanta = 0;
555 
556  // next section binomially assigns quanta to excitons and ions
557  G4int NumExcitons = BinomFluct(NumQuanta, fExcitationRatio / (1 + fExcitationRatio));
558  G4int NumIons = NumQuanta - NumExcitons;
559 
560  // this section calculates recombination following the modified Birks'
561  // Law of Doke, deposition by deposition, and may be overridden later
562  // in code if a low enough energy necessitates switching to the
563  // Thomas-Imel box model for recombination instead (determined by site)
564  G4double dE, dx = 0, LET = 0, recombProb;
565  dE = TotalEnergyDeposit / CLHEP::MeV;
566  if (particleName != "e-" && particleName != "e+" && z1 != z2 && particleName != "mu-" &&
567  particleName != "mu+") {
568  //in other words, if it's a gamma,ion,proton,alpha,pion,et al. do not
569  //use the step length provided by Geant4 because it's not relevant,
570  //instead calculate an estimated LET and range of the electrons that
571  //would have been produced if Geant4 could track them
572  LET = CalculateElectronLET(1000 * dE, z1);
573  if (LET) dx = dE / (Density * LET); //find the range based on the LET
574  if (abs(aParticle->GetPDGcode()) == 2112) dx = 0;
575  }
576  else { //normal case of an e-/+ energy deposition recorded by Geant
577  dx = aStep.GetStepLength() / CLHEP::cm;
578  if (dx) LET = (dE / dx) * (1 / Density); //lin. energy xfer (prop. to dE/dx)
579  if (LET > 0 && dE > 0 && dx > 0) {
580  G4double ratio = CalculateElectronLET(dE * 1e3, z1) / LET;
581  if (j == 1 && ratio < 0.7 && !ThomasImelTail && particleName == "e-") {
582  dx /= ratio;
583  LET *= ratio;
584  }
585  }
586  }
587  DokeBirks[1] = DokeBirks[0] / (1 - DokeBirks[2]); //B=A/(1-C) (see paper)
588  //Doke/Birks' Law as spelled out in the NEST paper
589  recombProb = (DokeBirks[0] * LET) / (1 + DokeBirks[1] * LET) + DokeBirks[2];
590  if (Phase == kStateLiquid) {
591  if (z1 == 54) recombProb *= (Density / Density_LXe);
592  if (z1 == 18) recombProb *= (Density / Density_LAr);
593  }
594  //check against unphysicality resulting from rounding errors
595  if (recombProb < 0) recombProb = 0;
596  if (recombProb > 1) recombProb = 1;
597  //use binomial distribution to assign photons, electrons, where photons
598  //are excitons plus recombined ionization electrons, while final
599  //collected electrons are the "escape" (non-recombined) electrons
600  G4int NumPhotons = NumExcitons + BinomFluct(NumIons, recombProb);
601  G4int NumElectrons = NumQuanta - NumPhotons;
602 
603  fEnergyDep = TotalEnergyDeposit;
604  fNumIonElectrons = NumElectrons;
605  fNumScintPhotons = NumPhotons;
606 
607  // next section increments the numbers of excitons, ions, photons, and
608  // electrons for the appropriate interaction site; it only appears to
609  // be redundant by saving seemingly no longer needed exciton and ion
610  // counts, these having been already used to calculate the number of ph
611  // and e- above, whereas it does need this later for Thomas-Imel model
612  char numExc[80];
613  char numIon[80];
614  char numPho[80];
615  char numEle[80];
616  sprintf(numExc, "N_EXC_%i", counter);
617  sprintf(numIon, "N_ION_%i", counter);
618  aMaterialPropertiesTable->AddConstProperty(numExc, NumExcitons);
619  aMaterialPropertiesTable->AddConstProperty(numIon, NumIons);
620  sprintf(numPho, "N_PHO_%i", counter);
621  sprintf(numEle, "N_ELE_%i", counter);
622  aMaterialPropertiesTable->AddConstProperty(numPho, NumPhotons);
623  aMaterialPropertiesTable->AddConstProperty(numEle, NumElectrons);
624 
625  // increment and save the total track length, and save interaction
626  // times for later, when generating the scintillation quanta
627  char trackL[80];
628  char time00[80];
629  char time01[80];
630  char energy[80];
631  sprintf(trackL, "TRACK_%i", counter);
632  sprintf(energy, "ENRGY_%i", counter);
633  sprintf(time00, "TIME0_%i", counter);
634  sprintf(time01, "TIME1_%i", counter);
635  delta = aMaterialPropertiesTable->GetConstProperty(trackL);
636  G4double energ = aMaterialPropertiesTable->GetConstProperty(energy);
637  delta += dx * CLHEP::cm;
638  energ += dE * CLHEP::MeV;
639  aMaterialPropertiesTable->AddConstProperty(trackL, delta);
640  aMaterialPropertiesTable->AddConstProperty(energy, energ);
641  if (TotalEnergyDeposit > 0) {
642  G4double deltaTime = aMaterialPropertiesTable->GetConstProperty(time00);
643  //for charged particles, which continuously lose energy, use initial
644  //interaction time as the minimum time, otherwise use only the final
645  if (aParticle->GetCharge() != 0) {
646  if (t0 < deltaTime) aMaterialPropertiesTable->AddConstProperty(time00, t0);
647  }
648  else {
649  if (t1 < deltaTime) aMaterialPropertiesTable->AddConstProperty(time00, t1);
650  }
651  deltaTime = aMaterialPropertiesTable->GetConstProperty(time01);
652  //find the maximum possible scintillation "birth" time
653  if (t1 > deltaTime) aMaterialPropertiesTable->AddConstProperty(time01, t1);
654  }
655 
656  // begin the process of setting up creation of scint./ionization
657  TotalEnergyDeposit =
658  aMaterialPropertiesTable->GetConstProperty("ENERGY_DEPOSIT_TOT"); //get the total E deposited
659  InitialKinetEnergy =
660  aMaterialPropertiesTable->GetConstProperty("ENERGY_DEPOSIT_GOL"); //E that should have been
661  if (InitialKinetEnergy > HIENLIM && abs(aParticle->GetPDGcode()) != 2112) fVeryHighEnergy = true;
662  G4double safety; //margin of error for TotalE.. - InitialKinetEnergy
663  if (fVeryHighEnergy && !fExcitedNucleus)
664  safety = 0.2 * CLHEP::keV;
665  else
666  safety = 2. * CLHEP::eV;
667 
668  //force a scintillation dump for NR and for full nuclear de-excitation
669  if (!anExcitationEnergy && pDef->GetParticleType() == "nucleus" &&
670  aTrack.GetTrackStatus() != fAlive && !fAlpha)
671  InitialKinetEnergy = TotalEnergyDeposit;
672  if (particleName == "neutron" || particleName == "antineutron")
673  InitialKinetEnergy = TotalEnergyDeposit;
674 
675  //force a dump of all saved scintillation under the following
676  //conditions: energy goal reached, and current particle dead, or an
677  //error has occurred and total has exceeded goal (shouldn't happen)
678  if (std::abs(TotalEnergyDeposit - InitialKinetEnergy) < safety ||
679  TotalEnergyDeposit >= InitialKinetEnergy) {
680  dx = 0;
681  dE = 0;
682  //calculate the total number of quanta from all sites and all
683  //interactions so that the number of secondaries gets set correctly
684  NumPhotons = 0;
685  NumElectrons = 0;
686  for (i = 0; i < j; i++) {
687  sprintf(numPho, "N_PHO_%d", i);
688  sprintf(numEle, "N_ELE_%d", i);
689  sprintf(trackL, "TRACK_%d", i);
690  sprintf(energy, "ENRGY_%d", i);
691  //add up track lengths of all sites, for a total LET calc (later)
692  dx += aMaterialPropertiesTable->GetConstProperty(trackL);
693  dE += aMaterialPropertiesTable->GetConstProperty(energy);
694  }
695 
696  G4int buffer = 100;
697  if (fVeryHighEnergy) buffer = 1;
698  fParticleChange.SetNumberOfSecondaries(buffer * (NumPhotons + NumElectrons));
699 
700  if (fTrackSecondariesFirst) {
701  if (aTrack.GetTrackStatus() == fAlive) fParticleChange.ProposeTrackStatus(fSuspend);
702  }
703 
704  // begin the loop over all sites which generates all the quanta
705  for (i = 0; i < j; i++) {
706  // get the position X,Y,Z, exciton and ion numbers, total track
707  // length of the site, and interaction times
708  sprintf(xCoord, "POS_X_%d", i);
709  sprintf(yCoord, "POS_Y_%d", i);
710  sprintf(zCoord, "POS_Z_%d", i);
711  sprintf(numExc, "N_EXC_%d", i);
712  sprintf(numIon, "N_ION_%d", i);
713  sprintf(numPho, "N_PHO_%d", i);
714  sprintf(numEle, "N_ELE_%d", i);
715  NumExcitons = (G4int)aMaterialPropertiesTable->GetConstProperty(numExc);
716  NumIons = (G4int)aMaterialPropertiesTable->GetConstProperty(numIon);
717  sprintf(trackL, "TRACK_%d", i);
718  sprintf(energy, "ENRGY_%d", i);
719  sprintf(time00, "TIME0_%d", i);
720  sprintf(time01, "TIME1_%d", i);
721  delta = aMaterialPropertiesTable->GetConstProperty(trackL);
722  energ = aMaterialPropertiesTable->GetConstProperty(energy);
723  t0 = aMaterialPropertiesTable->GetConstProperty(time00);
724  t1 = aMaterialPropertiesTable->GetConstProperty(time01);
725 
726  //if site is small enough, override the Doke/Birks' model with
727  //Thomas-Imel, but not if we're dealing with super-high energy
728  //particles, and if it's NR force Thomas-Imel (though NR should be
729  //already short enough in track even up to O(100) keV)
730  if ((delta < R0 && !fVeryHighEnergy) || z2 == z1 || fAlpha) {
731  if (z1 == 54 && ElectricField && //see NEST paper for ER formula
732  Phase == kStateLiquid) {
733  if (abs(z1 - z2) && //electron recoil
734  abs(aParticle->GetPDGcode()) != 2112) {
735  ThomasImel = 0.056776 * pow(ElectricField, -0.11844);
736  if (fAlpha) //technically ER, but special
737  ThomasImel = 0.057675 * pow(ElectricField, -0.49362);
738  } //end electron recoil (ER)
739  else { //nuclear recoil
740  // spline of NR data of C.E. Dahl PhD Thesis Princeton '09
741  // functions found using zunzun.com
742  ThomasImel = -0.15169 * pow((ElectricField + 215.25), 0.01811) + 0.20952;
743  } //end NR information
744  // Never let LY exceed 0-field yield!
745  if (ThomasImel > 0.19) ThomasImel = 0.19;
746  if (ThomasImel < 0.00) ThomasImel = 0.00;
747  } //end non-zero E-field segment
748  if (Phase == kStateLiquid) {
749  if (z1 == 54) ThomasImel *= pow((Density / 2.84), 0.3);
750  if (z1 == 18) ThomasImel *= pow((Density / Density_LAr), 0.3);
751  }
752  //calculate the Thomas-Imel recombination probability, which
753  //depends on energy via NumIons, but not on dE/dx, and protect
754  //against seg fault by ensuring a positive number of ions
755  if (NumIons > 0) {
756  G4double xi;
757  xi = (G4double(NumIons) / 4.) * ThomasImel;
758  if (InitialKinetEnergy == 9.4 * CLHEP::keV) {
759  G4double NumIonsEff =
760  1.066e7 * pow(t0 / CLHEP::ns, -1.303) * (0.17163 + 162.32 / (ElectricField + 191.39));
761  if (NumIonsEff > 1e6) NumIonsEff = 1e6;
762  xi = (G4double(NumIonsEff) / 4.) * ThomasImel;
763  }
764  if (fKr83m && ElectricField == 0) xi = (G4double(1.3 * NumIons) / 4.) * ThomasImel;
765  recombProb = 1 - log(1 + xi) / xi;
766  if (recombProb < 0) recombProb = 0;
767  if (recombProb > 1) recombProb = 1;
768  }
769  //just like Doke: simple binomial distribution
770  NumPhotons = NumExcitons + BinomFluct(NumIons, recombProb);
771  NumElectrons = (NumExcitons + NumIons) - NumPhotons;
772  //override Doke NumPhotons and NumElectrons
773  aMaterialPropertiesTable->AddConstProperty(numPho, NumPhotons);
774  aMaterialPropertiesTable->AddConstProperty(numEle, NumElectrons);
775  }
776 
777  // grab NumPhotons/NumElectrons, which come from Birks if
778  // the Thomas-Imel block of code above was not executed
779  NumPhotons = (G4int)aMaterialPropertiesTable->GetConstProperty(numPho);
780  NumElectrons = (G4int)aMaterialPropertiesTable->GetConstProperty(numEle);
781  // extra Fano factor caused by recomb. fluct.
782  G4double FanoFactor = 0; //ionization channel
783  if (Phase == kStateLiquid && fYieldFactor == 1) {
784  FanoFactor = 2575.9 * pow((ElectricField + 15.154), -0.64064) - 1.4707;
785  if (fKr83m) TotalEnergyDeposit = 4 * CLHEP::keV;
786  if ((dE / CLHEP::keV) <= 100 && ElectricField >= 0) {
787  G4double keVee = (TotalEnergyDeposit / (100. * CLHEP::keV));
788  if (keVee <= 0.06)
789  FanoFactor *= -0.00075 + 0.4625 * keVee + 34.375 * pow(keVee, 2.);
790  else
791  FanoFactor *= 0.069554 + 1.7322 * keVee - .80215 * pow(keVee, 2.);
792  }
793  }
794  if (Phase == kStateGas && Density > 0.5)
795  FanoFactor = 0.42857 - 4.7857 * Density + 7.8571 * pow(Density, 2.);
796  if (FanoFactor <= 0 || fVeryHighEnergy) FanoFactor = 0;
797  NumQuanta = NumPhotons + NumElectrons;
798  if (z1 == 54 && FanoFactor)
799  NumElectrons =
800  G4int(floor(GaussGen.fire(NumElectrons, sqrt(FanoFactor * NumElectrons)) + 0.5));
801  NumPhotons = NumQuanta - NumElectrons;
802  if (NumElectrons <= 0) NumElectrons = 0;
803  if (NumPhotons <= 0)
804  NumPhotons = 0;
805  else { //other effects
806  if (fAlpha) //bi-excitonic quenching due to high dE/dx
807  NumPhotons = BinomFluct(NumPhotons, biExc);
808  NumPhotons = BinomFluct(NumPhotons, QE_EFF);
809  }
810  NumElectrons = G4int(floor(NumElectrons * phe_per_e + 0.5));
811 
812  // new stuff to make Kr-83m work properly
813  if (fKr83m || InitialKinetEnergy == 9.4 * CLHEP::keV) fKr83m += dE / CLHEP::keV;
814  if (fKr83m > 41) fKr83m = 0;
815  if (SinglePhase) //for a 1-phase det. don't propagate e-'s
816  NumElectrons = 0; //saves simulation time
817 
818  // reset material properties numExc, numIon, numPho, numEle, as
819  // their values have been used or stored elsewhere already
820  aMaterialPropertiesTable->AddConstProperty(numExc, 0);
821  aMaterialPropertiesTable->AddConstProperty(numIon, 0);
822  aMaterialPropertiesTable->AddConstProperty(numPho, 0);
823  aMaterialPropertiesTable->AddConstProperty(numEle, 0);
824 
825  // start particle creation loop
826  if (InitialKinetEnergy < MAX_ENE && InitialKinetEnergy > MIN_ENE && !fMultipleScattering)
827  NumQuanta = NumPhotons + NumElectrons;
828  else
829  NumQuanta = 0;
830  for (k = 0; k < NumQuanta; k++) {
831  G4double sampledEnergy;
832  std::unique_ptr<G4DynamicParticle> aQuantum;
833 
834  // Generate random direction
835  G4double cost = 1. - 2. * UniformGen.fire();
836  G4double sint = std::sqrt((1. - cost) * (1. + cost));
837  G4double phi = CLHEP::twopi * UniformGen.fire();
838  G4double sinp = std::sin(phi);
839  G4double cosp = std::cos(phi);
840  G4double px = sint * cosp;
841  G4double py = sint * sinp;
842  G4double pz = cost;
843 
844  // Create momentum direction vector
845  G4ParticleMomentum photonMomentum(px, py, pz);
846 
847  // case of photon-specific stuff
848  if (k < NumPhotons) {
849  // Determine polarization of new photon
850  G4double sx = cost * cosp;
851  G4double sy = cost * sinp;
852  G4double sz = -sint;
853  G4ThreeVector photonPolarization(sx, sy, sz);
854  G4ThreeVector perp = photonMomentum.cross(photonPolarization);
855  phi = CLHEP::twopi * UniformGen.fire();
856  sinp = std::sin(phi);
857  cosp = std::cos(phi);
858  photonPolarization = cosp * photonPolarization + sinp * perp;
859  photonPolarization = photonPolarization.unit();
860 
861  // Generate a new photon or electron:
862  sampledEnergy = GaussGen.fire(PhotMean, PhotWidth);
863  aQuantum =
864  std::make_unique<G4DynamicParticle>(G4OpticalPhoton::OpticalPhoton(), photonMomentum);
865  aQuantum->SetPolarization(
866  photonPolarization.x(), photonPolarization.y(), photonPolarization.z());
867  }
868 
869  else { // this else statement is for ionization electrons
870  if (ElectricField) {
871  // point all electrons straight up, for drifting
872  G4ParticleMomentum electronMomentum(0, 0, -FieldSign);
873  aQuantum = std::make_unique<G4DynamicParticle>(G4ThermalElectron::ThermalElectron(),
874  electronMomentum);
875  if (Phase == kStateGas) {
876  sampledEnergy = GetGasElectronDriftSpeed(ElectricField, nDensity);
877  }
878  else
879  sampledEnergy =
880  GetLiquidElectronDriftSpeed(Temperature, ElectricField, MillerDriftSpeed, z1);
881  }
882  else {
883  // use "photonMomentum" for the electrons in the case of zero
884  // electric field, which is just randomized vector we made
885  aQuantum = std::make_unique<G4DynamicParticle>(G4ThermalElectron::ThermalElectron(),
886  photonMomentum);
887  sampledEnergy = 1.38e-23 * (CLHEP::joule / CLHEP::kelvin) * Temperature;
888  }
889  }
890 
891  //assign energy to make particle real
892  aQuantum->SetKineticEnergy(sampledEnergy);
893 
894  // Generate new G4Track object:
895  // emission time distribution
896 
897  // first an initial birth time is provided that is typically
898  // <<1 ns after the initial interaction in the simulation, then
899  // singlet, triplet lifetimes, and recombination time, are
900  // handled here, to create a realistic S1 pulse shape/timing
901  G4double aSecondaryTime = t0 + UniformGen.fire() * (t1 - t0) + evtStrt;
902  if (tau1 < 0) { tau1 = 0; }
903  if (tau3 < 0) { tau3 = 0; }
904  if (tauR < 0) { tauR = 0; }
905  if (aQuantum->GetDefinition()->GetParticleName() == "opticalphoton") {
906  if (abs(z2 - z1) && !fAlpha && //electron recoil
907  abs(aParticle->GetPDGcode()) != 2112) {
908  LET = (energ / CLHEP::MeV) / (delta / CLHEP::cm) * (1 / Density); //avg LET over all
909  //in future, this will be done interaction by interaction
910  // Next, find the recombination time, which is LET-dependent
911  // via ionization density (Kubota et al. Phys. Rev. B 20
912  // (1979) 3486). We find the LET-dependence by fitting to the
913  // E-dependence (Akimov et al. Phys. Lett. B 524 (2002) 245).
914  if (Phase == kStateLiquid && z1 == 54)
915  tauR =
916  3.5 * ((1 + 0.41 * LET) / (0.18 * LET)) * CLHEP::ns * exp(-0.00900 * ElectricField);
917  //field dependence based on fitting Fig. 9 of Dawson et al.
918  //NIM A 545 (2005) 690
919  //singlet-triplet ratios adapted from Kubota 1979, converted
920  //into correct units to work here, and separately done for
921  //excitation and recombination processes for electron recoils
922  //and assumed same for all LET (may vary)
923  SingTripRatioX = GaussGen.fire(0.17, 0.05);
924  SingTripRatioR = GaussGen.fire(0.8, 0.2);
925  if (z1 == 18) {
926  SingTripRatioR =
927  0.2701 + 0.003379 * LET - 4.7338e-5 * pow(LET, 2.) + 8.1449e-6 * pow(LET, 3.);
928  SingTripRatioX = SingTripRatioR;
929  if (LET < 3) {
930  SingTripRatioX = GaussGen.fire(0.36, 0.06);
931  SingTripRatioR = GaussGen.fire(0.5, 0.2);
932  }
933  }
934  }
935  else if (fAlpha) { //alpha particles
936  SingTripRatioR = GaussGen.fire(2.3, 0.51);
937  //currently based on Dawson 05 and Tey. 11 (arXiv:1103.3689)
938  //real ratio is likely a gentle function of LET
939  if (z1 == 18)
940  SingTripRatioR = (-0.065492 + 1.9996 * exp(-dE / CLHEP::MeV)) /
941  (1 + 0.082154 / pow(dE / CLHEP::MeV, 2.)) +
942  2.1811;
943  SingTripRatioX = SingTripRatioR;
944  }
945  else { //nuclear recoil
946  //based loosely on Hitachi et al. Phys. Rev. B 27 (1983) 5279
947  //with an eye to reproducing Akimov 2002 Fig. 9
948  SingTripRatioR = GaussGen.fire(7.8, 1.5);
949  if (z1 == 18) SingTripRatioR = 0.22218 * pow(energ / CLHEP::keV, 0.48211);
950  SingTripRatioX = SingTripRatioR;
951  }
952  // now, use binomial distributions to determine singlet and
953  // triplet states (and do separately for initially excited guys
954  // and recombining)
955  if (k > NumExcitons) {
956  //the recombination time is non-exponential, but approximates
957  //to exp at long timescales (see Kubota '79)
958  aSecondaryTime += tauR * (1. / UniformGen.fire() - 1);
959  if (UniformGen.fire() < SingTripRatioR / (1 + SingTripRatioR))
960  aSecondaryTime -= tau1 * log(UniformGen.fire());
961  else
962  aSecondaryTime -= tau3 * log(UniformGen.fire());
963  }
964  else {
965  if (UniformGen.fire() < SingTripRatioX / (1 + SingTripRatioX))
966  aSecondaryTime -= tau1 * log(UniformGen.fire());
967  else
968  aSecondaryTime -= tau3 * log(UniformGen.fire());
969  }
970  }
971  else { //electron trapping at the liquid/gas interface
972  G4double gainField = 12;
973  G4double tauTrap = 884.83 - 62.069 * gainField;
974  if (Phase == kStateLiquid) aSecondaryTime -= tauTrap * CLHEP::ns * log(UniformGen.fire());
975  }
976 
977  // emission position distribution --
978  // Generate the position of a new photon or electron, with NO
979  // stochastic variation because that could lead to particles
980  // being mistakenly generated outside of your active region by
981  // Geant4, but real-life finite detector position resolution
982  // wipes out any effects from here anyway...
983  std::string sx(xCoord);
984  std::string sy(yCoord);
985  std::string sz(zCoord);
986  x0[0] = aMaterialPropertiesTable->GetConstProperty(xCoord);
987  x0[1] = aMaterialPropertiesTable->GetConstProperty(yCoord);
988  x0[2] = aMaterialPropertiesTable->GetConstProperty(zCoord);
989  G4double radius = sqrt(pow(x0[0], 2.) + pow(x0[1], 2.));
990  //re-scale radius to ensure no generation of quanta outside
991  //the active volume of your simulation due to Geant4 rounding
992  if (radius >= R_TOL) {
993  if (x0[0] == 0) { x0[0] = 1 * CLHEP::nm; }
994  if (x0[1] == 0) { x0[1] = 1 * CLHEP::nm; }
995  radius -= R_TOL;
996  phi = atan(x0[1] / x0[0]);
997  x0[0] = fabs(radius * cos(phi)) * ((fabs(x0[0])) / (x0[0]));
998  x0[1] = fabs(radius * sin(phi)) * ((fabs(x0[1])) / (x0[1]));
999  }
1000  //position of the new secondary particle is ready for use
1001  G4ThreeVector aSecondaryPosition = x0;
1002  if (k >= NumPhotons && diffusion && ElectricField > 0) {
1003  G4double D_T = 64 * pow(1e-3 * ElectricField, -.17);
1004  //fit to Aprile and Doke 2009, arXiv:0910.4956 (Fig. 12)
1005  G4double D_L = 13.859 * pow(1e-3 * ElectricField, -0.58559);
1006  //fit to Aprile and Doke and Sorensen 2011, arXiv:1102.2865
1007  if (Phase == kStateLiquid && z1 == 18) {
1008  D_T = 93.342 * pow(ElectricField / nDensity, 0.041322);
1009  D_L = 0.15 * D_T;
1010  }
1011  if (Phase == kStateGas && z1 == 54) {
1012  D_L = 4.265 + 19097 / ElectricField - 1.7397e6 / pow(ElectricField, 2.) +
1013  1.2477e8 / pow(ElectricField, 3.);
1014  D_T *= 0.01;
1015  }
1016  D_T *= CLHEP::cm2 / CLHEP::s;
1017  D_L *= CLHEP::cm2 / CLHEP::s;
1018  if (ElectricField < 100 && Phase == kStateLiquid) D_L = 8 * CLHEP::cm2 / CLHEP::s;
1019  G4double vDrift = sqrt((2 * sampledEnergy) / (EMASS));
1020  if (BORDER == 0) x0[2] = 0;
1021  G4double sigmaDT = sqrt(2 * D_T * fabs(BORDER - x0[2]) / vDrift);
1022  G4double sigmaDL = sqrt(2 * D_L * fabs(BORDER - x0[2]) / vDrift);
1023  G4double dr = std::abs(GaussGen.fire(0., sigmaDT));
1024  phi = CLHEP::twopi * UniformGen.fire();
1025  aSecondaryPosition[0] += cos(phi) * dr;
1026  aSecondaryPosition[1] += sin(phi) * dr;
1027  aSecondaryPosition[2] += GaussGen.fire(0., sigmaDL);
1028  radius =
1029  std::sqrt(std::pow(aSecondaryPosition[0], 2.) + std::pow(aSecondaryPosition[1], 2.));
1030  if (aSecondaryPosition[2] >= BORDER && Phase == kStateLiquid) {
1031  if (BORDER != 0) aSecondaryPosition[2] = BORDER - R_TOL;
1032  }
1033  if (aSecondaryPosition[2] <= PMT && !GlobalFields) aSecondaryPosition[2] = PMT + R_TOL;
1034  } //end of electron diffusion code
1035 
1036  // GEANT4 business: stuff you need to make a new track
1037  if (aSecondaryTime < 0) aSecondaryTime = 0; //no neg. time
1038  /*
1039  auto aSecondaryTrack =
1040  std::make_unique<G4Track>(aQuantum,aSecondaryTime,aSecondaryPosition);
1041  if ( k < NumPhotons || radius < R_MAX )
1042  {
1043 
1044  */
1045  }
1046 
1047  //reset bunch of things when done with an interaction site
1048  aMaterialPropertiesTable->AddConstProperty(xCoord, 999 * CLHEP::km);
1049  aMaterialPropertiesTable->AddConstProperty(yCoord, 999 * CLHEP::km);
1050  aMaterialPropertiesTable->AddConstProperty(zCoord, 999 * CLHEP::km);
1051  aMaterialPropertiesTable->AddConstProperty(trackL, 0 * CLHEP::um);
1052  aMaterialPropertiesTable->AddConstProperty(energy, 0 * CLHEP::eV);
1053  aMaterialPropertiesTable->AddConstProperty(time00, DBL_MAX);
1054  aMaterialPropertiesTable->AddConstProperty(time01, -1 * CLHEP::ns);
1055 
1056  } //end of interaction site loop
1057 
1058  //more things to reset...
1059  aMaterialPropertiesTable->AddConstProperty("TOTALNUM_INT_SITES", 0);
1060  aMaterialPropertiesTable->AddConstProperty("ENERGY_DEPOSIT_TOT", 0 * CLHEP::keV);
1061  aMaterialPropertiesTable->AddConstProperty("ENERGY_DEPOSIT_GOL", 0 * CLHEP::MeV);
1062  fExcitedNucleus = false;
1063  fAlpha = false;
1064  }
1065 
1066  //don't do anything when you're not ready to scintillate
1067  else {
1068  fParticleChange.SetNumberOfSecondaries(0);
1069  return fParticleChange;
1070  }
1071 
1072  return fParticleChange;
1073 }
1074 
1075 //----------------------------------------------------------------------------
1076 G4double NestAlg::GetGasElectronDriftSpeed(G4double /*efieldinput*/, G4double /*density*/)
1077 {
1078  std::cout << "WARNING: NestAlg::GetGasElectronDriftSpeed(G4double, G4double) "
1079  << "is not defined, returning bogus value of -999." << std::endl;
1080 
1081  return -999.;
1082 }
1083 
1084 //----------------------------------------------------------------------------
1085 G4double NestAlg::GetLiquidElectronDriftSpeed(G4double tempinput,
1086  G4double efieldinput,
1087  G4bool Miller,
1088  G4int Z)
1089 {
1090  if (efieldinput < 0) efieldinput *= (-1);
1091  //Liquid equation one (165K) coefficients
1092  G4double onea = 144623.235704015, oneb = 850.812714257629, onec = 1192.87056676815,
1093  oned = -395969.575204061, onef = -355.484170008875, oneg = -227.266219627672,
1094  oneh = 223831.601257495, onei = 6.1778950907965, onej = 18.7831533426398,
1095  onek = -76132.6018884368;
1096  //Liquid equation two (200K) coefficients
1097  G4double twoa = 17486639.7118995, twob = -113.174284723134, twoc = 28.005913193763,
1098  twod = 167994210.094027, twof = -6766.42962575088, twog = 901.474643115395,
1099  twoh = -185240292.471665, twoi = -633.297790813084, twoj = 87.1756135457949;
1100  //Liquid equation three (230K) coefficients
1101  G4double thra = 10626463726.9833, thrb = 224025158.134792, thrc = 123254826.300172,
1102  thrd = -4563.5678061122, thrf = -1715.269592063, thrg = -694181.921834368,
1103  thrh = -50.9753281079838, thri = 58.3785811395493, thrj = 201512.080026704;
1104  G4double y1 = 0, y2 = 0, f1 = 0, f2 = 0, f3 = 0, edrift = 0, t1 = 0, t2 = 0, slope = 0,
1105  intercept = 0;
1106 
1107  //Equations defined
1108  f1 = onea / (1 + exp(-(efieldinput - oneb) / onec)) +
1109  oned / (1 + exp(-(efieldinput - onef) / oneg)) +
1110  oneh / (1 + exp(-(efieldinput - onei) / onej)) + onek;
1111  f2 = twoa / (1 + exp(-(efieldinput - twob) / twoc)) +
1112  twod / (1 + exp(-(efieldinput - twof) / twog)) +
1113  twoh / (1 + exp(-(efieldinput - twoi) / twoj));
1114  f3 = thra * exp(-thrb * efieldinput) + thrc * exp(-(pow(efieldinput - thrd, 2)) / (thrf * thrf)) +
1115  thrg * exp(-(pow(efieldinput - thrh, 2) / (thri * thri))) + thrj;
1116 
1117  if (efieldinput < 20 && efieldinput >= 0) {
1118  f1 = 2951 * efieldinput;
1119  f2 = 5312 * efieldinput;
1120  f3 = 7101 * efieldinput;
1121  }
1122  //Cases for tempinput decides which 2 equations to use lin. interpolation
1123  if (tempinput < 200.0 && tempinput > 165.0) {
1124  y1 = f1;
1125  y2 = f2;
1126  t1 = 165.0;
1127  t2 = 200.0;
1128  }
1129  if (tempinput < 230.0 && tempinput > 200.0) {
1130  y1 = f2;
1131  y2 = f3;
1132  t1 = 200.0;
1133  t2 = 230.0;
1134  }
1135  if ((tempinput > 230.0 || tempinput < 165.0) && !Miller) {
1136  G4cout << "\nWARNING: TEMPERATURE OUT OF RANGE (165-230 K)\n";
1137  return 0;
1138  }
1139  if (tempinput == 165.0)
1140  edrift = f1;
1141  else if (tempinput == 200.0)
1142  edrift = f2;
1143  else if (tempinput == 230.0)
1144  edrift = f3;
1145  else { //Linear interpolation
1146  //frac=(tempinput-t1)/(t2-t1);
1147  slope = (y1 - y2) / (t1 - t2);
1148  intercept = y1 - slope * t1;
1149  edrift = slope * tempinput + intercept;
1150  }
1151 
1152  if (Miller) {
1153  if (efieldinput <= 40.)
1154  edrift = -0.13274 + 0.041082 * efieldinput - 0.0006886 * pow(efieldinput, 2.) +
1155  5.5503e-6 * pow(efieldinput, 3.);
1156  else
1157  edrift = 0.060774 * efieldinput / pow(1 + 0.11336 * pow(efieldinput, 0.5218), 2.);
1158  if (efieldinput >= 1e5) edrift = 2.7;
1159  if (efieldinput >= 100)
1160  edrift -= 0.017 * (tempinput - 163);
1161  else
1162  edrift += 0.017 * (tempinput - 163);
1163  edrift *= 1e5; //put into units of cm/sec. from mm/usec.
1164  }
1165  if (Z == 18)
1166  edrift = 1e5 * (.097384 * pow(log10(efieldinput), 3.0622) - .018614 * sqrt(efieldinput));
1167  if (edrift < 0) edrift = 0.;
1168  edrift = 0.5 * EMASS * pow(edrift * CLHEP::cm / CLHEP::s, 2.);
1169  return edrift;
1170 }
1171 
1172 //----------------------------------------------------------------------------
1173 G4double NestAlg::CalculateElectronLET(G4double E, G4int Z)
1174 {
1175  G4double LET;
1176  switch (Z) {
1177  case 54:
1178  //use a spline fit to online ESTAR data
1179  if (E >= 1)
1180  LET = 58.482 - 61.183 * log10(E) + 19.749 * pow(log10(E), 2) + 2.3101 * pow(log10(E), 3) -
1181  3.3469 * pow(log10(E), 4) + 0.96788 * pow(log10(E), 5) - 0.12619 * pow(log10(E), 6) +
1182  0.0065108 * pow(log10(E), 7);
1183  //at energies <1 keV, use a different spline, determined manually by
1184  //generating sub-keV electrons in Geant4 and looking at their ranges, since
1185  //ESTAR does not go this low
1186  else if (E > 0 && E < 1)
1187  LET = 6.9463 + 815.98 * E - 4828 * pow(E, 2) + 17079 * pow(E, 3) - 36394 * pow(E, 4) +
1188  44553 * pow(E, 5) - 28659 * pow(E, 6) + 7483.8 * pow(E, 7);
1189  else
1190  LET = 0;
1191  break;
1192  case 18:
1193  default:
1194  if (E >= 1)
1195  LET = 116.70 - 162.97 * log10(E) + 99.361 * pow(log10(E), 2) - 33.405 * pow(log10(E), 3) +
1196  6.5069 * pow(log10(E), 4) - 0.69334 * pow(log10(E), 5) + .031563 * pow(log10(E), 6);
1197  else if (E > 0 && E < 1)
1198  LET = 100;
1199  else
1200  LET = 0;
1201  }
1202  return LET;
1203 }
1204 
1205 //----------------------------------------------------------------------------
1206 G4int NestAlg::BinomFluct(G4int N0, G4double prob)
1207 {
1208  CLHEP::RandGauss GaussGen(fEngine);
1209  CLHEP::RandFlat UniformGen(fEngine);
1210 
1211  G4double mean = N0 * prob;
1212  G4double sigma = sqrt(N0 * prob * (1 - prob));
1213  G4int N1 = 0;
1214  if (prob == 0.00) return N1;
1215  if (prob == 1.00) return N0;
1216 
1217  if (N0 < 10) {
1218  for (G4int i = 0; i < N0; i++) {
1219  if (UniformGen.fire() < prob) N1++;
1220  }
1221  }
1222  else {
1223  N1 = G4int(floor(GaussGen.fire(mean, sigma) + 0.5));
1224  }
1225  if (N1 > N0) N1 = N0;
1226  if (N1 < 0) N1 = 0;
1227  return N1;
1228 }
1229 
1230 //----------------------------------------------------------------------------
1231 void NestAlg::InitMatPropValues(G4MaterialPropertiesTable* nobleElementMat, int z)
1232 {
1233  char xCoord[80];
1234  char yCoord[80];
1235  char zCoord[80];
1236  char numExc[80];
1237  char numIon[80];
1238  char numPho[80];
1239  char numEle[80];
1240  char trackL[80];
1241  char time00[80];
1242  char time01[80];
1243  char energy[80];
1244 
1245  // for loop to initialize the interaction site mat'l properties
1246  for (G4int i = 0; i < 10000; i++) {
1247  sprintf(xCoord, "POS_X_%d", i);
1248  sprintf(yCoord, "POS_Y_%d", i);
1249  sprintf(zCoord, "POS_Z_%d", i);
1250  nobleElementMat->AddConstProperty(xCoord, 999 * CLHEP::km);
1251  nobleElementMat->AddConstProperty(yCoord, 999 * CLHEP::km);
1252  nobleElementMat->AddConstProperty(zCoord, 999 * CLHEP::km);
1253  sprintf(numExc, "N_EXC_%d", i);
1254  sprintf(numIon, "N_ION_%d", i);
1255  sprintf(numPho, "N_PHO_%d", i);
1256  sprintf(numEle, "N_ELE_%d", i);
1257  nobleElementMat->AddConstProperty(numExc, 0);
1258  nobleElementMat->AddConstProperty(numIon, 0);
1259  nobleElementMat->AddConstProperty(numPho, 0);
1260  nobleElementMat->AddConstProperty(numEle, 0);
1261  sprintf(trackL, "TRACK_%d", i);
1262  sprintf(energy, "ENRGY_%d", i);
1263  sprintf(time00, "TIME0_%d", i);
1264  sprintf(time01, "TIME1_%d", i);
1265  nobleElementMat->AddConstProperty(trackL, 0 * CLHEP::um);
1266  nobleElementMat->AddConstProperty(energy, 0 * CLHEP::eV);
1267  nobleElementMat->AddConstProperty(time00, DBL_MAX);
1268  nobleElementMat->AddConstProperty(time01, -1 * CLHEP::ns);
1269  }
1270 
1271  // we initialize the total number of interaction sites, a variable for
1272  // updating the amount of energy deposited thus far in the medium, and a
1273  // variable for storing the amount of energy expected to be deposited
1274  nobleElementMat->AddConstProperty("TOTALNUM_INT_SITES", 0);
1275  nobleElementMat->AddConstProperty("ENERGY_DEPOSIT_TOT", 0 * CLHEP::keV);
1276  nobleElementMat->AddConstProperty("ENERGY_DEPOSIT_GOL", 0 * CLHEP::MeV);
1277 
1278  fElementPropInit[z] = true;
1279 
1280  return;
1281 }
1282 
1283 //----------------------------------------------------------------------------
1284 G4double UnivScreenFunc(G4double E, G4double Z, G4double A)
1285 {
1286  G4double a_0 = 5.29e-11 * CLHEP::m;
1287  G4double a = 0.626 * a_0 * pow(Z, (-1. / 3.));
1288  G4double epsilon_0 = 8.854e-12 * (CLHEP::farad / CLHEP::m);
1289  G4double epsilon =
1290  (a * E * 2. * CLHEP::twopi * epsilon_0) / (2 * pow(CLHEP::eplus, 2.) * pow(Z, 2.));
1291  G4double zeta_0 = pow(Z, (1. / 6.));
1292  G4double m_N = A * 1.66e-27 * CLHEP::kg;
1293  G4double hbar = 6.582e-16 * CLHEP::eV * CLHEP::s;
1294  if (Z == 54) {
1295  epsilon *= 1.068; //zeta_0 = 1.63;
1296  } //special case for LXe from Bezrukov et al. 2011
1297  G4double s_n = log(1 + 1.1383 * epsilon) /
1298  (2. * (epsilon + 0.01321 * pow(epsilon, 0.21226) + 0.19593 * sqrt(epsilon)));
1299  G4double s_e =
1300  (a_0 * zeta_0 / a) * hbar *
1301  sqrt(8 * epsilon * 2. * CLHEP::twopi * epsilon_0 / (a * m_N * pow(CLHEP::eplus, 2.)));
1302  return 1.38e5 * 0.5 * (1 + tanh(50 * epsilon - 0.25)) * epsilon * (s_e / s_n);
1303 }
code to link reconstructed objects back to the MC truth information
#define AVO
Definition: NestAlg.cxx:1
#define Density_LAr
Definition: NestAlg.cxx:26
G4bool ThomasImelTail
Definition: NestAlg.cxx:42
TTree * t1
Definition: plottest35.C:26
Float_t y1[n_points_granero]
Definition: compare.C:5
#define R_TOL
Definition: NestAlg.cxx:23
#define EMASS
Definition: NestAlg.cxx:2
#define BORDER
Definition: NestAlg.cxx:6
G4double GetLiquidElectronDriftSpeed(double T, double F, G4bool M, G4int Z)
Definition: NestAlg.cxx:1085
Float_t x1[n_points_granero]
Definition: compare.C:5
G4double biExc
Definition: NestAlg.cxx:44
NestAlg(CLHEP::HepRandomEngine &engine)
Definition: NestAlg.cxx:49
G4bool SinglePhase
Definition: NestAlg.cxx:42
Double_t z
Definition: plot.C:276
#define HIENLIM
Definition: NestAlg.cxx:22
#define WIN
Definition: NestAlg.cxx:12
constexpr auto abs(T v)
Returns the absolute value of the argument.
double fYieldFactor
turns scint. on/off
Definition: NestAlg.h:39
Float_t f2
G4VParticleChange fParticleChange
pointer to G4VParticleChange
Definition: NestAlg.h:45
int fNumIonElectrons
number of ionization electrons produced by step
Definition: NestAlg.h:43
kilovolt_as<> kilovolt
Type of potential stored in kilovolt, in double precision.
G4int BinomFluct(G4int N0, G4double prob)
Definition: NestAlg.cxx:1206
Float_t y2[n_points_geant4]
Definition: compare.C:26
#define PMT
Definition: NestAlg.cxx:19
Float_t Z
Definition: plot.C:37
int fNumScintPhotons
number of photons produced by the step
Definition: NestAlg.h:42
#define a2
Float_t f3
auto counter(T begin, T end)
Returns an object to iterate values from begin to end in a range-for loop.
Definition: counter.h:295
static G4ThermalElectron * ThermalElectron()
Float_t E
Definition: plot.C:20
Float_t radius
Definition: plot.C:23
G4bool diffusion
Definition: NestAlg.cxx:40
const G4VParticleChange & CalculateIonizationAndScintillation(G4Track const &aTrack, G4Step const &aStep)
Definition: NestAlg.cxx:81
G4double UnivScreenFunc(G4double E, G4double Z, G4double A)
double energy
Definition: plottest35.C:25
#define Density_LXe
Definition: NestAlg.cxx:25
Float_t f1
G4double GetGasElectronDriftSpeed(G4double efieldinput, G4double density)
Definition: NestAlg.cxx:1076
CLHEP::HepRandomEngine & fEngine
random engine
Definition: NestAlg.h:49
G4double CalculateElectronLET(G4double E, G4int Z)
Definition: NestAlg.cxx:1173
TTree * t2
Definition: plottest35.C:36
#define GAT
Definition: NestAlg.cxx:16
#define phe_per_e
Definition: NestAlg.cxx:9
#define SRF
Definition: NestAlg.cxx:15
double mean(const std::vector< short > &wf, size_t start, size_t nsample)
Definition: UtilFunc.cxx:13
#define QE_EFF
Definition: NestAlg.cxx:8
void InitMatPropValues(G4MaterialPropertiesTable *nobleElementMat, int z)
Definition: NestAlg.cxx:1231
#define BOT
Definition: NestAlg.cxx:18
#define MillerDriftSpeed
Definition: NestAlg.cxx:3
#define ANE
Definition: NestAlg.cxx:14
double fExcitationRatio
Definition: NestAlg.h:40
#define CTH
Definition: NestAlg.cxx:17
Float_t e
Definition: plot.C:35
G4bool OutElectrons
Definition: NestAlg.cxx:42
#define TOP
Definition: NestAlg.cxx:13
#define MIN_ENE
Definition: NestAlg.cxx:20
std::map< int, bool > fElementPropInit
Definition: NestAlg.h:46
double fEnergyDep
energy deposited by the step
Definition: NestAlg.h:44
#define a1