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