LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
Public Member Functions | Private Member Functions | Private Attributes | List of all members
NestAlg Class Reference

#include "NestAlg.h"

Public Member Functions

 NestAlg (CLHEP::HepRandomEngine &engine)
 
 NestAlg (double yieldFactor, CLHEP::HepRandomEngine &engine)
 
const G4VParticleChange & CalculateIonizationAndScintillation (G4Track const &aTrack, G4Step const &aStep)
 
void SetScintillationYieldFactor (double const &yf)
 
void SetScintillationExcitationRatio (double const &er)
 
int NumberScintillationPhotons () const
 
int NumberIonizationElectrons () const
 
double EnergyDeposition () const
 

Private Member Functions

G4double GetGasElectronDriftSpeed (G4double efieldinput, G4double density)
 
G4double GetLiquidElectronDriftSpeed (double T, double F, G4bool M, G4int Z)
 
G4double CalculateElectronLET (G4double E, G4int Z)
 
G4double UnivScreenFunc (G4double E, G4double Z, G4double A)
 
G4int BinomFluct (G4int N0, G4double prob)
 
void InitMatPropValues (G4MaterialPropertiesTable *nobleElementMat, int z)
 

Private Attributes

double fYieldFactor
 turns scint. on/off More...
 
double fExcitationRatio
 
int fNumScintPhotons
 number of photons produced by the step More...
 
int fNumIonElectrons
 number of ionization electrons produced by step More...
 
double fEnergyDep
 energy deposited by the step More...
 
G4VParticleChange fParticleChange
 pointer to G4VParticleChange More...
 
std::map< int, bool > fElementPropInit
 
CLHEP::HepRandomEngine & fEngine
 random engine More...
 

Detailed Description

Definition at line 16 of file NestAlg.h.

Constructor & Destructor Documentation

NestAlg::NestAlg ( CLHEP::HepRandomEngine &  engine)

Definition at line 49 of file NestAlg.cxx.

References fElementPropInit.

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 }
NestAlg::fYieldFactor
double fYieldFactor
turns scint. on/off
Definition: NestAlg.h:39
NestAlg::fNumIonElectrons
int fNumIonElectrons
number of ionization electrons produced by step
Definition: NestAlg.h:43
NestAlg::fNumScintPhotons
int fNumScintPhotons
number of photons produced by the step
Definition: NestAlg.h:42
NestAlg::fEngine
CLHEP::HepRandomEngine & fEngine
random engine
Definition: NestAlg.h:49
NestAlg::fExcitationRatio
double fExcitationRatio
Definition: NestAlg.h:40
NestAlg::fElementPropInit
std::map< int, bool > fElementPropInit
Definition: NestAlg.h:46
NestAlg::fEnergyDep
double fEnergyDep
energy deposited by the step
Definition: NestAlg.h:44
NestAlg::NestAlg ( double  yieldFactor,
CLHEP::HepRandomEngine &  engine 
)

Definition at line 65 of file NestAlg.cxx.

References fElementPropInit.

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 }
NestAlg::fYieldFactor
double fYieldFactor
turns scint. on/off
Definition: NestAlg.h:39
NestAlg::fNumIonElectrons
int fNumIonElectrons
number of ionization electrons produced by step
Definition: NestAlg.h:43
NestAlg::fNumScintPhotons
int fNumScintPhotons
number of photons produced by the step
Definition: NestAlg.h:42
NestAlg::fEngine
CLHEP::HepRandomEngine & fEngine
random engine
Definition: NestAlg.h:49
NestAlg::fExcitationRatio
double fExcitationRatio
Definition: NestAlg.h:40
NestAlg::fElementPropInit
std::map< int, bool > fElementPropInit
Definition: NestAlg.h:46
NestAlg::fEnergyDep
double fEnergyDep
energy deposited by the step
Definition: NestAlg.h:44

Member Function Documentation

G4int NestAlg::BinomFluct ( G4int  N0,
G4double  prob 
)
private

Definition at line 1206 of file NestAlg.cxx.

References fEngine, and pmtana::mean().

Referenced by CalculateIonizationAndScintillation().

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 }
NestAlg::fEngine
CLHEP::HepRandomEngine & fEngine
random engine
Definition: NestAlg.h:49
pmtana::mean
double mean(const std::vector< short > &wf, size_t start, size_t nsample)
Definition: UtilFunc.cxx:13
G4double NestAlg::CalculateElectronLET ( G4double  E,
G4int  Z 
)
private

Definition at line 1173 of file NestAlg.cxx.

Referenced by CalculateIonizationAndScintillation().

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 }
Z
Float_t Z
Definition: plot.C:37
E
Float_t E
Definition: plot.C:20
const G4VParticleChange & NestAlg::CalculateIonizationAndScintillation ( G4Track const &  aTrack,
G4Step const &  aStep 
)

Definition at line 81 of file NestAlg.cxx.

References a1, a2, util::abs(), ANE, AVO, biExc, BinomFluct(), BORDER, BOT, CalculateElectronLET(), util::counter(), CTH, Density_LAr, Density_LXe, diffusion, e, EMASS, energy, fElementPropInit, fEnergyDep, fEngine, fExcitationRatio, fNumIonElectrons, fNumScintPhotons, fParticleChange, fYieldFactor, GAT, GetGasElectronDriftSpeed(), GetLiquidElectronDriftSpeed(), HIENLIM, InitMatPropValues(), MillerDriftSpeed, MIN_ENE, OutElectrons, phe_per_e, PMT, QE_EFF, R_TOL, radius, SinglePhase, SRF, t1, G4ThermalElectron::ThermalElectron(), ThomasImelTail, TOP, WIN, and x1.

Referenced by larg4::ISCalculationNEST::CalculateIonizationAndScintillation().

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 }
t0
code to link reconstructed objects back to the MC truth information
Definition: DirectHitParticleAssns_tool.cc:18
AVO
#define AVO
Definition: NestAlg.cxx:1
Density_LAr
#define Density_LAr
Definition: NestAlg.cxx:26
ThomasImelTail
G4bool ThomasImelTail
Definition: NestAlg.cxx:42
t1
TTree * t1
Definition: plottest35.C:26
R_TOL
#define R_TOL
Definition: NestAlg.cxx:23
EMASS
#define EMASS
Definition: NestAlg.cxx:2
BORDER
#define BORDER
Definition: NestAlg.cxx:6
NestAlg::GetLiquidElectronDriftSpeed
G4double GetLiquidElectronDriftSpeed(double T, double F, G4bool M, G4int Z)
Definition: NestAlg.cxx:1085
x1
Float_t x1[n_points_granero]
Definition: compare.C:5
biExc
G4double biExc
Definition: NestAlg.cxx:44
SinglePhase
G4bool SinglePhase
Definition: NestAlg.cxx:42
HIENLIM
#define HIENLIM
Definition: NestAlg.cxx:22
WIN
#define WIN
Definition: NestAlg.cxx:12
util::abs
constexpr auto abs(T v)
Returns the absolute value of the argument.
Definition: constexpr_math.h:40
NestAlg::fYieldFactor
double fYieldFactor
turns scint. on/off
Definition: NestAlg.h:39
NestAlg::fParticleChange
G4VParticleChange fParticleChange
pointer to G4VParticleChange
Definition: NestAlg.h:45
NestAlg::fNumIonElectrons
int fNumIonElectrons
number of ionization electrons produced by step
Definition: NestAlg.h:43
util::quantities::kilovolt
kilovolt_as<> kilovolt
Type of potential stored in kilovolt, in double precision.
Definition: electromagnetism.h:190
NestAlg::BinomFluct
G4int BinomFluct(G4int N0, G4double prob)
Definition: NestAlg.cxx:1206
PMT
#define PMT
Definition: NestAlg.cxx:19
NestAlg::fNumScintPhotons
int fNumScintPhotons
number of photons produced by the step
Definition: NestAlg.h:42
a2
#define a2
Definition: HoughBaseAlg.cxx:64
util::counter
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
G4ThermalElectron::ThermalElectron
static G4ThermalElectron * ThermalElectron()
Definition: G4ThermalElectron.cxx:102
radius
Float_t radius
Definition: plot.C:23
diffusion
G4bool diffusion
Definition: NestAlg.cxx:40
energy
double energy
Definition: plottest35.C:25
Density_LXe
#define Density_LXe
Definition: NestAlg.cxx:25
NestAlg::GetGasElectronDriftSpeed
G4double GetGasElectronDriftSpeed(G4double efieldinput, G4double density)
Definition: NestAlg.cxx:1076
NestAlg::fEngine
CLHEP::HepRandomEngine & fEngine
random engine
Definition: NestAlg.h:49
NestAlg::CalculateElectronLET
G4double CalculateElectronLET(G4double E, G4int Z)
Definition: NestAlg.cxx:1173
GAT
#define GAT
Definition: NestAlg.cxx:16
phe_per_e
#define phe_per_e
Definition: NestAlg.cxx:9
SRF
#define SRF
Definition: NestAlg.cxx:15
QE_EFF
#define QE_EFF
Definition: NestAlg.cxx:8
NestAlg::InitMatPropValues
void InitMatPropValues(G4MaterialPropertiesTable *nobleElementMat, int z)
Definition: NestAlg.cxx:1231
BOT
#define BOT
Definition: NestAlg.cxx:18
MillerDriftSpeed
#define MillerDriftSpeed
Definition: NestAlg.cxx:3
ANE
#define ANE
Definition: NestAlg.cxx:14
NestAlg::fExcitationRatio
double fExcitationRatio
Definition: NestAlg.h:40
CTH
#define CTH
Definition: NestAlg.cxx:17
e
Float_t e
Definition: plot.C:35
OutElectrons
G4bool OutElectrons
Definition: NestAlg.cxx:42
TOP
#define TOP
Definition: NestAlg.cxx:13
MIN_ENE
#define MIN_ENE
Definition: NestAlg.cxx:20
NestAlg::fElementPropInit
std::map< int, bool > fElementPropInit
Definition: NestAlg.h:46
NestAlg::fEnergyDep
double fEnergyDep
energy deposited by the step
Definition: NestAlg.h:44
a1
#define a1
Definition: HoughBaseAlg.cxx:63
double NestAlg::EnergyDeposition ( ) const
inline

Definition at line 29 of file NestAlg.h.

References E, UnivScreenFunc(), Z, and z.

Referenced by larg4::ISCalculationNEST::CalculateIonizationAndScintillation().

29 { return fEnergyDep; }
NestAlg::fEnergyDep
double fEnergyDep
energy deposited by the step
Definition: NestAlg.h:44
G4double NestAlg::GetGasElectronDriftSpeed ( G4double  efieldinput,
G4double  density 
)
private

Definition at line 1076 of file NestAlg.cxx.

Referenced by CalculateIonizationAndScintillation().

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 }
G4double NestAlg::GetLiquidElectronDriftSpeed ( double  T,
double  F,
G4bool  M,
G4int  Z 
)
private

Definition at line 1085 of file NestAlg.cxx.

References EMASS, f1, f2, f3, t1, t2, y1, and y2.

Referenced by CalculateIonizationAndScintillation().

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 }
t1
TTree * t1
Definition: plottest35.C:26
y1
Float_t y1[n_points_granero]
Definition: compare.C:5
EMASS
#define EMASS
Definition: NestAlg.cxx:2
f2
Float_t f2
Definition: comparison_ascii.C:56
y2
Float_t y2[n_points_geant4]
Definition: compare.C:26
Z
Float_t Z
Definition: plot.C:37
f3
Float_t f3
Definition: comparison_ascii.C:56
f1
Float_t f1
Definition: comparison_ascii.C:56
t2
TTree * t2
Definition: plottest35.C:36
void NestAlg::InitMatPropValues ( G4MaterialPropertiesTable *  nobleElementMat,
int  z 
)
private

Definition at line 1231 of file NestAlg.cxx.

References energy, fElementPropInit, and z.

Referenced by CalculateIonizationAndScintillation().

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 }
z
Double_t z
Definition: plot.C:276
energy
double energy
Definition: plottest35.C:25
NestAlg::fElementPropInit
std::map< int, bool > fElementPropInit
Definition: NestAlg.h:46
int NestAlg::NumberIonizationElectrons ( ) const
inline

Definition at line 28 of file NestAlg.h.

Referenced by larg4::ISCalculationNEST::CalculateIonizationAndScintillation().

28 { return fNumIonElectrons; }
NestAlg::fNumIonElectrons
int fNumIonElectrons
number of ionization electrons produced by step
Definition: NestAlg.h:43
int NestAlg::NumberScintillationPhotons ( ) const
inline

Definition at line 27 of file NestAlg.h.

Referenced by larg4::ISCalculationNEST::CalculateIonizationAndScintillation().

27 { return fNumScintPhotons; }
NestAlg::fNumScintPhotons
int fNumScintPhotons
number of photons produced by the step
Definition: NestAlg.h:42
void NestAlg::SetScintillationExcitationRatio ( double const &  er)
inline

Definition at line 26 of file NestAlg.h.

26 { fExcitationRatio = er; }
NestAlg::fExcitationRatio
double fExcitationRatio
Definition: NestAlg.h:40
void NestAlg::SetScintillationYieldFactor ( double const &  yf)
inline

Definition at line 25 of file NestAlg.h.

25 { fYieldFactor = yf; }
NestAlg::fYieldFactor
double fYieldFactor
turns scint. on/off
Definition: NestAlg.h:39
G4double NestAlg::UnivScreenFunc ( G4double  E,
G4double  Z,
G4double  A 
)
private

Member Data Documentation

std::map<int, bool> NestAlg::fElementPropInit
private

map of noble element z to flag for whether that element's material properties table has been initialized

Definition at line 46 of file NestAlg.h.

Referenced by CalculateIonizationAndScintillation(), InitMatPropValues(), and NestAlg().

double NestAlg::fEnergyDep
private

energy deposited by the step

Definition at line 44 of file NestAlg.h.

Referenced by CalculateIonizationAndScintillation().

CLHEP::HepRandomEngine& NestAlg::fEngine
private

random engine

Definition at line 49 of file NestAlg.h.

Referenced by BinomFluct(), and CalculateIonizationAndScintillation().

double NestAlg::fExcitationRatio
private

N_ex/N_i, the dimensionless ratio of initial excitons to ions

Definition at line 40 of file NestAlg.h.

Referenced by CalculateIonizationAndScintillation().

int NestAlg::fNumIonElectrons
private

number of ionization electrons produced by step

Definition at line 43 of file NestAlg.h.

Referenced by CalculateIonizationAndScintillation().

int NestAlg::fNumScintPhotons
private

number of photons produced by the step

Definition at line 42 of file NestAlg.h.

Referenced by CalculateIonizationAndScintillation().

G4VParticleChange NestAlg::fParticleChange
private

pointer to G4VParticleChange

Definition at line 45 of file NestAlg.h.

Referenced by CalculateIonizationAndScintillation().

double NestAlg::fYieldFactor
private

turns scint. on/off

Definition at line 39 of file NestAlg.h.

Referenced by CalculateIonizationAndScintillation().

The documentation for this class was generated from the following files: