LArSoft  v07_13_02
Liquid Argon Software toolkit - http://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 12 of file NestAlg.h.

Constructor & Destructor Documentation

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

Definition at line 50 of file NestAlg.cxx.

References fElementPropInit.

51  : fYieldFactor(0)
52  , fExcitationRatio(0)
53  , fNumScintPhotons(0)
54  , fNumIonElectrons(0)
55  , fEnergyDep(0.)
56  , fEngine(engine)
57 {
58  fElementPropInit[2] = false;
59  fElementPropInit[10] = false;
60  fElementPropInit[18] = false;
61  fElementPropInit[36] = false;
62  fElementPropInit[54] = false;
63 }
NestAlg::fYieldFactor
double fYieldFactor
turns scint. on/off
Definition: NestAlg.h:37
NestAlg::fNumIonElectrons
int fNumIonElectrons
number of ionization electrons produced by step
Definition: NestAlg.h:41
NestAlg::fNumScintPhotons
int fNumScintPhotons
number of photons produced by the step
Definition: NestAlg.h:40
NestAlg::fEngine
CLHEP::HepRandomEngine & fEngine
random engine
Definition: NestAlg.h:47
NestAlg::fElementPropInit
std::map< int, bool > fElementPropInit
Definition: NestAlg.h:44
NestAlg::fExcitationRatio
double fExcitationRatio
Definition: NestAlg.h:38
NestAlg::fEnergyDep
double fEnergyDep
energy deposited by the step
Definition: NestAlg.h:42
NestAlg::NestAlg ( double  yieldFactor,
CLHEP::HepRandomEngine &  engine 
)

Definition at line 66 of file NestAlg.cxx.

References fElementPropInit.

67  : fYieldFactor(yieldFactor)
68  , fExcitationRatio(0.)
69  , fNumScintPhotons(0)
70  , fNumIonElectrons(0)
71  , fEnergyDep(0.)
72  , fEngine(engine)
73 {
74  fElementPropInit[2] = false;
75  fElementPropInit[10] = false;
76  fElementPropInit[18] = false;
77  fElementPropInit[36] = false;
78  fElementPropInit[54] = false;
79 }
NestAlg::fYieldFactor
double fYieldFactor
turns scint. on/off
Definition: NestAlg.h:37
NestAlg::fNumIonElectrons
int fNumIonElectrons
number of ionization electrons produced by step
Definition: NestAlg.h:41
NestAlg::fNumScintPhotons
int fNumScintPhotons
number of photons produced by the step
Definition: NestAlg.h:40
NestAlg::fEngine
CLHEP::HepRandomEngine & fEngine
random engine
Definition: NestAlg.h:47
NestAlg::fElementPropInit
std::map< int, bool > fElementPropInit
Definition: NestAlg.h:44
NestAlg::fExcitationRatio
double fExcitationRatio
Definition: NestAlg.h:38
NestAlg::fEnergyDep
double fEnergyDep
energy deposited by the step
Definition: NestAlg.h:42

Member Function Documentation

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

Definition at line 1190 of file NestAlg.cxx.

References fEngine, and pmtana::mean().

Referenced by CalculateIonizationAndScintillation().

1190  {
1191  CLHEP::RandGauss GaussGen(fEngine);
1192  CLHEP::RandFlat UniformGen(fEngine);
1193 
1194  G4double mean = N0*prob;
1195  G4double sigma = sqrt(N0*prob*(1-prob));
1196  G4int N1 = 0;
1197  if ( prob == 0.00 ) return N1;
1198  if ( prob == 1.00 ) return N0;
1199 
1200  if ( N0 < 10 ) {
1201  for(G4int i = 0; i < N0; i++) {
1202  if(UniformGen.fire() < prob) N1++;
1203  }
1204  }
1205  else {
1206  N1 = G4int(floor(GaussGen.fire(mean,sigma)+0.5));
1207  }
1208  if ( N1 > N0 ) N1 = N0;
1209  if ( N1 < 0 ) N1 = 0;
1210  return N1;
1211 }
NestAlg::fEngine
CLHEP::HepRandomEngine & fEngine
random engine
Definition: NestAlg.h:47
pmtana::mean
double mean(const std::vector< short > &wf, size_t start, size_t nsample)
Definition: UtilFunc.cxx:15
G4double NestAlg::CalculateElectronLET ( G4double  E,
G4int  Z 
)
private

Definition at line 1162 of file NestAlg.cxx.

Referenced by CalculateIonizationAndScintillation().

1162  {
1163  G4double LET;
1164  switch ( Z ) {
1165  case 54:
1166  //use a spline fit to online ESTAR data
1167  if ( E >= 1 ) LET = 58.482-61.183*log10(E)+19.749*pow(log10(E),2)+
1168  2.3101*pow(log10(E),3)-3.3469*pow(log10(E),4)+
1169  0.96788*pow(log10(E),5)-0.12619*pow(log10(E),6)+0.0065108*pow(log10(E),7);
1170  //at energies <1 keV, use a different spline, determined manually by
1171  //generating sub-keV electrons in Geant4 and looking at their ranges, since
1172  //ESTAR does not go this low
1173  else if ( E>0 && E<1 ) LET = 6.9463+815.98*E-4828*pow(E,2)+17079*pow(E,3)-
1174  36394*pow(E,4)+44553*pow(E,5)-28659*pow(E,6)+7483.8*pow(E,7);
1175  else
1176  LET = 0;
1177  break;
1178  case 18: default:
1179  if ( E >= 1 ) LET = 116.70-162.97*log10(E)+99.361*pow(log10(E),2)-
1180  33.405*pow(log10(E),3)+6.5069*pow(log10(E),4)-
1181  0.69334*pow(log10(E),5)+.031563*pow(log10(E),6);
1182  else if ( E>0 && E<1 ) LET = 100;
1183  else
1184  LET = 0;
1185  }
1186  return LET;
1187 }
E
Float_t E
Definition: plot.C:23
Z
Float_t Z
Definition: plot.C:39
const G4VParticleChange & NestAlg::CalculateIonizationAndScintillation ( G4Track const &  aTrack,
G4Step const &  aStep 
)

Definition at line 82 of file NestAlg.cxx.

References a1, a2, ANE, AVO, biExc, BinomFluct(), BORDER, BOT, CalculateElectronLET(), 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, s, SinglePhase, SRF, t1, G4ThermalElectron::ThermalElectron(), ThomasImelTail, TOP, WIN, and x1.

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

84 {
85  CLHEP::RandGauss GaussGen(fEngine);
86  CLHEP::RandFlat UniformGen(fEngine);
87 
88 
89  // reset the variables accessed by other objects
90  // make the energy deposit the energy in this step,
91  // set the number of electrons and photons to 0
92  fEnergyDep = aStep.GetTotalEnergyDeposit();
93  fNumIonElectrons = 0.;
94  fNumScintPhotons = 0.;
95 
96 
97  if ( !fYieldFactor ) //set YF=0 when you want S1Light off in your sim
98  return fParticleChange;
99 
100  bool fTrackSecondariesFirst = false;
101  bool fExcitedNucleus = false;
102  bool fVeryHighEnergy = false;
103  bool fAlpha = false;
104  bool fMultipleScattering = false;
105  double fKr83m = 0.;
106 
107  if( aTrack.GetParentID() == 0 && aTrack.GetCurrentStepNumber() == 1 ) {
108  fExcitedNucleus = false; //an initialization or reset
109  fVeryHighEnergy = false; //initializes or (later) resets this
110  fAlpha = false; //ditto
111  fMultipleScattering = false;
112  }
113 
114  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
115  G4ParticleDefinition *pDef = aParticle->GetDefinition();
116  G4String particleName = pDef->GetParticleName();
117  const G4Material* aMaterial = aStep.GetPreStepPoint()->GetMaterial();
118  const G4Material* bMaterial = aStep.GetPostStepPoint()->GetMaterial();
119 
120  if((particleName == "neutron" || particleName == "antineutron") &&
121  aStep.GetTotalEnergyDeposit() <= 0)
122  return fParticleChange;
123 
124  // code for determining whether the present/next material is noble
125  // element, or, in other words, for checking if either is a valid NEST
126  // scintillating material, and save Z for later L calculation, or
127  // return if no valid scintillators are found on this step, which is
128  // protection against G4Exception or seg. fault/violation
129  G4Element *ElementA = NULL, *ElementB = NULL;
130  if (aMaterial) {
131  const G4ElementVector* theElementVector1 =
132  aMaterial->GetElementVector();
133  ElementA = (*theElementVector1)[0];
134  }
135  if (bMaterial) {
136  const G4ElementVector* theElementVector2 =
137  bMaterial->GetElementVector();
138  ElementB = (*theElementVector2)[0];
139  }
140  G4int z1,z2,j=1; G4bool NobleNow=false,NobleLater=false;
141  if (ElementA) z1 = (G4int)(ElementA->GetZ()); else z1 = -1;
142  if (ElementB) z2 = (G4int)(ElementB->GetZ()); else z2 = -1;
143  if ( z1==2 || z1==10 || z1==18 || z1==36 || z1==54 ) {
144  NobleNow = true;
145  // j = (G4int)aMaterial->GetMaterialPropertiesTable()->
146  // GetConstProperty("TOTALNUM_INT_SITES"); //get current number
147  //if ( j < 0 ) {
148  if ( aTrack.GetParentID() == 0 && !fElementPropInit[z1] ) {
149  InitMatPropValues(aMaterial->GetMaterialPropertiesTable(), z1);
150  j = 0; //no sites yet
151  } //material properties initialized
152  } //end of atomic number check
153  if ( z2==2 || z2==10 || z2==18 || z2==36 || z2==54 ) {
154  NobleLater = true;
155  // j = (G4int)bMaterial->GetMaterialPropertiesTable()->
156  // GetConstProperty("TOTALNUM_INT_SITES");
157  //if ( j < 0 ) {
158  if ( aTrack.GetParentID() == 0 && !fElementPropInit[z2] ) {
159  InitMatPropValues(bMaterial->GetMaterialPropertiesTable(), z2);
160  j = 0; //no sites yet
161  } //material properties initialized
162  } //end of atomic number check
163 
164  if ( !NobleNow && !NobleLater )
165  return fParticleChange;
166 
167  // retrieval of the particle's position, time, attributes at both the
168  // beginning and the end of the current step along its track
169  G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
170  G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
171  G4ThreeVector x1 = pPostStepPoint->GetPosition();
172  G4ThreeVector x0 = pPreStepPoint->GetPosition();
173  G4double evtStrt = pPreStepPoint->GetGlobalTime();
174  G4double t0 = pPreStepPoint->GetLocalTime();
175  G4double t1 = pPostStepPoint->GetLocalTime();
176 
177  // now check if we're entering a scintillating material (inside) or
178  // leaving one (outside), in order to determine (later on in the code,
179  // based on the booleans inside & outside) whether to add/subtract
180  // energy that can potentially be deposited from the system
181  G4bool outside = false, inside = false, InsAndOuts = false;
182  G4MaterialPropertiesTable* aMaterialPropertiesTable =
183  aMaterial->GetMaterialPropertiesTable();
184  if ( NobleNow && !NobleLater ) outside = true;
185  if ( !NobleNow && NobleLater ) {
186  aMaterial = bMaterial; inside = true; z1 = z2;
187  aMaterialPropertiesTable = bMaterial->GetMaterialPropertiesTable();
188  }
189  if ( NobleNow && NobleLater &&
190  aMaterial->GetDensity() != bMaterial->GetDensity() )
191  InsAndOuts = true;
192 
193  // retrieve scintillation-related material properties
194  G4double Density = aMaterial->GetDensity()/(CLHEP::g/CLHEP::cm3);
195  G4double nDensity = Density*AVO; //molar mass factor applied below
196  G4int Phase = aMaterial->GetState(); //solid, liquid, or gas?
197  G4double ElectricField(0.), FieldSign(0.); //for field quenching of S1
198  G4bool GlobalFields = false;
199 
200  if ( (WIN == 0) && !TOP && !ANE && !SRF && !GAT && !CTH && !BOT && !PMT ) {
201  ElectricField = aMaterialPropertiesTable->GetConstProperty("ELECTRICFIELD");
202  GlobalFields = true;
203  }
204  else {
205  if ( x1[2] < WIN && x1[2] > TOP && Phase == kStateGas )
206  ElectricField = aMaterialPropertiesTable->GetConstProperty("ELECTRICFIELDWINDOW");
207  else if ( x1[2] < TOP && x1[2] > ANE && Phase == kStateGas )
208  ElectricField = aMaterialPropertiesTable->GetConstProperty("ELECTRICFIELDTOP");
209  else if ( x1[2] < ANE && x1[2] > SRF && Phase == kStateGas )
210  ElectricField = aMaterialPropertiesTable->
211  GetConstProperty("ELECTRICFIELDANODE");
212  else if ( x1[2] < SRF && x1[2] > GAT && Phase == kStateLiquid )
213  ElectricField = aMaterialPropertiesTable->
214  GetConstProperty("ELECTRICFIELDSURFACE");
215  else if ( x1[2] < GAT && x1[2] > CTH && Phase == kStateLiquid )
216  ElectricField = aMaterialPropertiesTable->
217  GetConstProperty("ELECTRICFIELDGATE");
218  else if ( x1[2] < CTH && x1[2] > BOT && Phase == kStateLiquid )
219  ElectricField = aMaterialPropertiesTable->
220  GetConstProperty("ELECTRICFIELDCATHODE");
221  else if ( x1[2] < BOT && x1[2] > PMT && Phase == kStateLiquid )
222  ElectricField = aMaterialPropertiesTable->
223  GetConstProperty("ELECTRICFIELDBOTTOM");
224  else
225  ElectricField = aMaterialPropertiesTable->
226  GetConstProperty("ELECTRICFIELD");
227  }
228  if ( ElectricField >= 0 ) FieldSign = 1; else FieldSign = -1;
229  ElectricField = std::abs((1e3*ElectricField)/(CLHEP::kilovolt/CLHEP::cm));
230  G4double Temperature = aMaterial->GetTemperature();
231  G4double ScintillationYield, ResolutionScale, R0 = 1.0*CLHEP::um,
232  DokeBirks[3], ThomasImel = 0.00, delta = 1*CLHEP::mm;
233  DokeBirks[0] = 0.00; DokeBirks[2] = 1.00;
234  G4double PhotMean = 7*CLHEP::eV, PhotWidth = 1.0*CLHEP::eV; //photon properties
235  G4double SingTripRatioR, SingTripRatioX, tau1, tau3, tauR = 0*CLHEP::ns;
236  switch ( z1 ) { //sort prop. by noble element atomic#
237  case 2: //helium
238  ScintillationYield = 1 / (41.3*CLHEP::eV); //all W's from noble gas book
239  fExcitationRatio = 0.00; //nominal (true value unknown)
240  ResolutionScale = 0.2; //Aprile, Bolotnikov, Bolozdynya, Doke
241  PhotMean = 15.9*CLHEP::eV;
242  tau1 = GaussGen.fire(10.0*CLHEP::ns,0.0*CLHEP::ns);
243  tau3 = 1.6e3*CLHEP::ns;
244  tauR = GaussGen.fire(13.00*CLHEP::s,2.00*CLHEP::s); //McKinsey et al. 2003
245  break;
246  case 10: //neon
247  ScintillationYield = 1 / (29.2*CLHEP::eV);
248  fExcitationRatio = 0.00; //nominal (true value unknown)
249  ResolutionScale = 0.13; //Aprile et. al book
250  PhotMean = 15.5*CLHEP::eV; 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  } nDensity /= 39.948; //molar mass in grams per mole
269  PhotMean = 9.69*CLHEP::eV; PhotWidth = 0.22*CLHEP::eV;
270  tau1 = GaussGen.fire(6.5*CLHEP::ns,0.8*CLHEP::ns); //err from wgted avg.
271  tau3 = GaussGen.fire(1300*CLHEP::ns,50*CLHEP::ns); //ibid.
272  tauR = GaussGen.fire(0.8*CLHEP::ns,0.2*CLHEP::ns); //Kubota 1979
273  biExc = 0.6; break;
274  case 36: //krypton
275  if ( Phase == kStateGas ) ScintillationYield = 1 / (30.0*CLHEP::eV);
276  else ScintillationYield = 1 / (15.0*CLHEP::eV);
277  fExcitationRatio = 0.08; //Aprile et. al book
278  ResolutionScale = 0.05; //Doke 1976
279  PhotMean = 8.43*CLHEP::eV;
280  tau1 = GaussGen.fire(2.8*CLHEP::ns,.04*CLHEP::ns);
281  tau3 = GaussGen.fire(93.*CLHEP::ns,1.1*CLHEP::ns);
282  tauR = GaussGen.fire(12.*CLHEP::ns,.76*CLHEP::ns);
283  break;
284  case 54: //xenon
285  default: nDensity /= 131.293;
286  ScintillationYield = 48.814+0.80877*Density+2.6721*pow(Density,2.);
287  ScintillationYield /= CLHEP::keV; //Units: [#quanta(ph/e-) per keV]
288  //W = 13.7 eV for all recoils, Dahl thesis (that's @2.84 g/cm^3)
289  //the exciton-to-ion ratio (~0.06 for LXe at 3 g/cm^3)
290  fExcitationRatio = 0.4-0.11131*Density-0.0026651*pow(Density,2.);
291  ResolutionScale = 1.00 * //Fano factor <<1
292  (0.12724-0.032152*Density-0.0013492*pow(Density,2.));
293  //~0.1 for GXe w/ formula from Bolotnikov et al. 1995
294  if ( Phase == kStateLiquid ) {
295  ResolutionScale *= 1.5; //to get it to be ~0.03 for LXe
296  R0 = 16.6*CLHEP::um; //for zero electric field
297  //length scale above which Doke model used instead of Thomas-Imel
298  if(ElectricField) //change it with field (see NEST paper)
299  R0 = 69.492*pow(ElectricField,-0.50422)*CLHEP::um;
300  if(ElectricField) { //formulae & values all from NEST paper
301  DokeBirks[0]= 19.171*pow(ElectricField+25.552,-0.83057)+0.026772;
302  DokeBirks[2] = 0.00; //only volume recombination (above)
303  }
304  else { //zero electric field magnitude
305  DokeBirks[0] = 0.18; //volume/columnar recombination factor (A)
306  DokeBirks[2] = 0.58; //geminate/Onsager recombination (C)
307  }
308  //"ThomasImel" is alpha/(a^2*v), the recomb. coeff.
309  ThomasImel = 0.05; //aka xi/(4*N_i) from the NEST paper
310  //distance used to determine when one is at a new interaction site
311  delta = 0.4*CLHEP::mm; //distance ~30 keV x-ray travels in LXe
312  PhotMean = 6.97*CLHEP::eV; PhotWidth = 0.23*CLHEP::eV;
313  // 178+/-14nmFWHM, taken from Jortner JchPh 42 '65.
314  //these singlet and triplet times may not be the ones you're
315  //used to, but are the world average: Kubota 79, Hitachi 83 (2
316  //data sets), Teymourian 11, Morikawa 89, and Akimov '02
317  tau1 = GaussGen.fire(3.1*CLHEP::ns,.7*CLHEP::ns); //err from wgted avg.
318  tau3 = GaussGen.fire(24.*CLHEP::ns,1.*CLHEP::ns); //ibid.
319  } //end liquid
320  else if ( Phase == kStateGas ) {
321  if(!fAlpha) fExcitationRatio=0.07; //Nygren NIM A 603 (2009) p. 340
322  else { biExc = 1.00;
323  ScintillationYield = 1 / (12.98*CLHEP::eV); } //Saito 2003
324  R0 = 0.0*CLHEP::um; //all Doke/Birks interactions (except for alphas)
325  G4double Townsend = (ElectricField/nDensity)*1e17;
326  DokeBirks[0] = 0.0000; //all geminate (except at zero, low fields)
327  DokeBirks[2] = 0.1933*pow(Density,2.6199)+0.29754 -
328  (0.045439*pow(Density,2.4689)+0.066034)*log10(ElectricField);
329  if ( ElectricField>6990 ) DokeBirks[2]=0.0;
330  if ( ElectricField<1000 ) DokeBirks[2]=0.2;
331  if ( ElectricField<100. ) { DokeBirks[0]=0.18; DokeBirks[2]=0.58; }
332  if( Density < 0.061 ) ThomasImel = 0.041973*pow(Density,1.8105);
333  else if( Density >= 0.061 && Density <= 0.167 )
334  ThomasImel=5.9583e-5+0.0048523*Density-0.023109*pow(Density,2.);
335  else ThomasImel = 6.2552e-6*pow(Density,-1.9963);
336  if(ElectricField)ThomasImel=1.2733e-5*pow(Townsend/Density,-0.68426);
337  // field\density dependence from Kobayashi 2004 and Saito 2003
338  PhotMean = 7.1*CLHEP::eV; PhotWidth = 0.2*CLHEP::eV;
339  tau1 = GaussGen.fire(5.18*CLHEP::ns,1.55*CLHEP::ns);
340  tau3 = GaussGen.fire(100.1*CLHEP::ns,7.9*CLHEP::ns);
341  } //end gas information (preliminary guesses)
342  else {
343  tau1 = 3.5*CLHEP::ns; tau3 = 20.*CLHEP::ns; tauR = 40.*CLHEP::ns;
344  } //solid Xe
345  }
346 
347  // log present and running tally of energy deposition in this section
348  G4double anExcitationEnergy = ((const G4Ions*)(pDef))->
349  GetExcitationEnergy(); //grab nuclear energy level
350  G4double TotalEnergyDeposit = //total energy deposited so far
351  aMaterialPropertiesTable->GetConstProperty( "ENERGY_DEPOSIT_TOT" );
352  G4bool convert = false, annihil = false;
353  //set up special cases for pair production and positron annihilation
354  if(pPreStepPoint->GetKineticEnergy()>=(2*CLHEP::electron_mass_c2) &&
355  !pPostStepPoint->GetKineticEnergy() &&
356  !aStep.GetTotalEnergyDeposit() && aParticle->GetPDGcode()==22) {
357  convert = true; TotalEnergyDeposit = CLHEP::electron_mass_c2;
358  }
359  if(pPreStepPoint->GetKineticEnergy() &&
360  !pPostStepPoint->GetKineticEnergy() &&
361  aParticle->GetPDGcode()==-11) {
362  annihil = true; TotalEnergyDeposit += aStep.GetTotalEnergyDeposit();
363  }
364  G4bool either = false;
365  if(inside || outside || convert || annihil || InsAndOuts) either=true;
366  //conditions for returning when energy deposits too low
367  if( anExcitationEnergy<100*CLHEP::eV && aStep.GetTotalEnergyDeposit()<1*CLHEP::eV &&
368  !either && !fExcitedNucleus )
369  return fParticleChange;
370  //add current deposit to total energy budget
371  if ( !annihil ) TotalEnergyDeposit += aStep.GetTotalEnergyDeposit();
372  if ( !convert ) aMaterialPropertiesTable->
373  AddConstProperty( "ENERGY_DEPOSIT_TOT", TotalEnergyDeposit );
374  //save current deposit for determining number of quanta produced now
375  TotalEnergyDeposit = aStep.GetTotalEnergyDeposit();
376 
377  // check what the current "goal" E is for dumping scintillation,
378  // often the initial kinetic energy of the parent particle, and deal
379  // with all other energy-related matters in this block of code
380  G4double InitialKinetEnergy = aMaterialPropertiesTable->
381  GetConstProperty( "ENERGY_DEPOSIT_GOL" );
382  //if zero, add up initial potential and kinetic energies now
383  if ( InitialKinetEnergy == 0 ) {
384  G4double tE = pPreStepPoint->GetKineticEnergy()+anExcitationEnergy;
385  if ( (fabs(tE-1.8*CLHEP::keV) < CLHEP::eV || fabs(tE-9.4*CLHEP::keV) < CLHEP::eV) &&
386  Phase == kStateLiquid && z1 == 54 ) tE = 9.4*CLHEP::keV;
387  if ( fKr83m && ElectricField != 0 )
388  DokeBirks[2] = 0.10;
389  aMaterialPropertiesTable->
390  AddConstProperty ( "ENERGY_DEPOSIT_GOL", tE );
391  //excited nucleus is special case where accuracy reduced for total
392  //energy deposition because of G4 inaccuracies and scintillation is
393  //forced-dumped when that nucleus is fully de-excited
394  if ( anExcitationEnergy ) fExcitedNucleus = true;
395  }
396  //if a particle is leaving, remove its kinetic energy from the goal
397  //energy, as this will never get deposited (if depositable)
398  if(outside){ aMaterialPropertiesTable->
399  AddConstProperty("ENERGY_DEPOSIT_GOL",
400  InitialKinetEnergy-pPostStepPoint->GetKineticEnergy());
401  if(aMaterialPropertiesTable->
402  GetConstProperty("ENERGY_DEPOSIT_GOL")<0)
403  aMaterialPropertiesTable->AddConstProperty("ENERGY_DEPOSIT_GOL",0);
404  }
405  //if a particle is coming back into your scintillator, then add its
406  //energy to the goal energy
407  if(inside) { aMaterialPropertiesTable->
408  AddConstProperty("ENERGY_DEPOSIT_GOL",
409  InitialKinetEnergy+pPreStepPoint->GetKineticEnergy());
410  if ( TotalEnergyDeposit > 0 && InitialKinetEnergy == 0 ) {
411  aMaterialPropertiesTable->AddConstProperty("ENERGY_DEPOSIT_GOL",0);
412  TotalEnergyDeposit = .000000;
413  }
414  }
415  if ( InsAndOuts ) {
416  //G4double dribble = pPostStepPoint->GetKineticEnergy() -
417  //pPreStepPoint->GetKineticEnergy();
418  aMaterialPropertiesTable->
419  AddConstProperty("ENERGY_DEPOSIT_GOL",(-0.1*CLHEP::keV)+
420  InitialKinetEnergy-pPostStepPoint->GetKineticEnergy());
421  InitialKinetEnergy = bMaterial->GetMaterialPropertiesTable()->
422  GetConstProperty("ENERGY_DEPOSIT_GOL");
423  bMaterial->GetMaterialPropertiesTable()->
424  AddConstProperty("ENERGY_DEPOSIT_GOL",(-0.1*CLHEP::keV)+
425  InitialKinetEnergy+pPreStepPoint->GetKineticEnergy());
426  if(aMaterialPropertiesTable->
427  GetConstProperty("ENERGY_DEPOSIT_GOL")<0)
428  aMaterialPropertiesTable->AddConstProperty("ENERGY_DEPOSIT_GOL",0);
429  if ( bMaterial->GetMaterialPropertiesTable()->
430  GetConstProperty("ENERGY_DEPOSIT_GOL") < 0 )
431  bMaterial->GetMaterialPropertiesTable()->
432  AddConstProperty ( "ENERGY_DEPOSIT_GOL", 0 );
433  }
434  InitialKinetEnergy = aMaterialPropertiesTable->
435  GetConstProperty("ENERGY_DEPOSIT_GOL"); //grab current goal E
436  if ( annihil ) //if an annihilation occurred, add energy of two gammas
437  InitialKinetEnergy += 2*CLHEP::electron_mass_c2;
438  //if pair production occurs, then subtract energy to cancel with the
439  //energy that will be added in the line above when the e+ dies
440  if ( convert )
441  InitialKinetEnergy -= 2*CLHEP::electron_mass_c2;
442  //update the relevant material property (goal energy)
443  aMaterialPropertiesTable->
444  AddConstProperty("ENERGY_DEPOSIT_GOL",InitialKinetEnergy);
445  if (anExcitationEnergy < 1e-100 && aStep.GetTotalEnergyDeposit()==0 &&
446  aMaterialPropertiesTable->GetConstProperty("ENERGY_DEPOSIT_GOL")==0 &&
447  aMaterialPropertiesTable->GetConstProperty("ENERGY_DEPOSIT_TOT")==0)
448  return fParticleChange;
449 
450  G4String procName;
451  if ( aTrack.GetCreatorProcess() )
452  procName = aTrack.GetCreatorProcess()->GetProcessName();
453  else
454  procName = "NULL";
455  if ( procName == "eBrem" && outside && !OutElectrons )
456  fMultipleScattering = true;
457 
458  // next 2 codeblocks deal with position-related things
459  if ( fAlpha ) delta = 1000.*CLHEP::km;
460  G4int i, k, counter = 0; G4double pos[3];
461  if ( outside ) { //leaving
462  if ( aParticle->GetPDGcode() == 11 && !OutElectrons )
463  fMultipleScattering = true;
464  x1 = x0; //prevents generation of quanta outside active volume
465  } //no scint. for e-'s that leave
466 
467  char xCoord[80]; char yCoord[80]; char zCoord[80];
468  G4bool exists = false; //for querying whether set-up of new site needed
469  for(i=0;i<j;i++) { //loop over all saved interaction sites
470  counter = i; //save site# for later use in storing properties
471  sprintf(xCoord,"POS_X_%d",i); sprintf(yCoord,"POS_Y_%d",i);
472  sprintf(zCoord,"POS_Z_%d",i);
473  pos[0] = aMaterialPropertiesTable->GetConstProperty(xCoord);
474  pos[1] = aMaterialPropertiesTable->GetConstProperty(yCoord);
475  pos[2] = aMaterialPropertiesTable->GetConstProperty(zCoord);
476  if ( sqrt(pow(x1[0]-pos[0],2.)+pow(x1[1]-pos[1],2.)+
477  pow(x1[2]-pos[2],2.)) < delta ) {
478  exists = true; break; //we find interaction is close to an old one
479  }
480  }
481  if(!exists && TotalEnergyDeposit) { //current interaction too far away
482  counter = j;
483  sprintf(xCoord,"POS_X_%i",j); sprintf(yCoord,"POS_Y_%i",j);
484  sprintf(zCoord,"POS_Z_%i",j);
485  //save 3-space coordinates of the new interaction site
486  aMaterialPropertiesTable->AddConstProperty( xCoord, x1[0] );
487  aMaterialPropertiesTable->AddConstProperty( yCoord, x1[1] );
488  aMaterialPropertiesTable->AddConstProperty( zCoord, x1[2] );
489  j++; //increment number of sites
490  aMaterialPropertiesTable-> //save
491  AddConstProperty( "TOTALNUM_INT_SITES", j );
492  }
493 
494  // this is where nuclear recoil "L" factor is handled: total yield is
495  // reduced for nuclear recoil as per Lindhard theory
496 
497  //we assume you have a mono-elemental scintillator only
498  //now, grab A's and Z's of current particle and of material (avg)
499  G4double a1 = ElementA->GetA();
500  z2 = pDef->GetAtomicNumber();
501  G4double a2 = (G4double)(pDef->GetAtomicMass());
502  if ( particleName == "alpha" || (z2 == 2 && a2 == 4) )
503  fAlpha = true; //used later to get S1 pulse shape correct for alpha
504  if ( fAlpha || abs(aParticle->GetPDGcode()) == 2112 )
505  a2 = a1; //get average A for element at hand
506  G4double epsilon = 11.5*(TotalEnergyDeposit/CLHEP::keV)*pow(z1,(-7./3.));
507  G4double gamma = 3.*pow(epsilon,0.15)+0.7*pow(epsilon,0.6)+epsilon;
508  G4double kappa = 0.133*pow(z1,(2./3.))*pow(a2,(-1./2.))*(2./3.);
509  //check if we are dealing with nuclear recoil (Z same as material)
510  if ( (z1 == z2 && pDef->GetParticleType() == "nucleus") ||
511  particleName == "neutron" || particleName == "antineutron" ) {
512  fYieldFactor=(kappa*gamma)/(1+kappa*gamma); //Lindhard factor
513  if ( z1 == 18 && Phase == kStateLiquid )
514  fYieldFactor=0.23*(1+exp(-5*epsilon)); //liquid argon L_eff
515  //just a few safety checks, like for recombProb below
516  if ( fYieldFactor > 1 ) fYieldFactor = 1;
517  if ( fYieldFactor < 0 ) fYieldFactor = 0;
518  if ( ElectricField == 0 && Phase == kStateLiquid ) {
519  if ( z1 == 54 ) ThomasImel = 0.19;
520  if ( z1 == 18 ) ThomasImel = 0.25;
521  } //special TIB parameters for nuclear recoil only, in LXe and LAr
522  fExcitationRatio = 0.69337 +
523  0.3065*exp(-0.008806*pow(ElectricField,0.76313));
524  }
525  else fYieldFactor = 1.000; //default
526 
527  // determine ultimate #quanta from current E-deposition (ph+e-)
528  G4double MeanNumberOfQuanta = //total mean number of exc/ions
529  ScintillationYield*TotalEnergyDeposit;
530  //the total number of either quanta produced is equal to product of the
531  //work function, the energy deposited, and yield reduction, for NR
532  G4double sigma = sqrt(ResolutionScale*MeanNumberOfQuanta); //Fano
533  G4int NumQuanta = //stochastic variation in NumQuanta
534  G4int(floor(GaussGen.fire(MeanNumberOfQuanta,sigma)+0.5));
535  G4double LeffVar = GaussGen.fire(fYieldFactor,0.25*fYieldFactor);
536  if (LeffVar > 1) { LeffVar = 1.00000; }
537  if (LeffVar < 0) { LeffVar = 0; }
538  if ( fYieldFactor < 1 ) NumQuanta = BinomFluct(NumQuanta,LeffVar);
539  //if E below work function, can't make any quanta, and if NumQuanta
540  //less than zero because Gaussian fluctuated low, update to zero
541  if(TotalEnergyDeposit < 1/ScintillationYield || NumQuanta < 0)
542  NumQuanta = 0;
543 
544  // next section binomially assigns quanta to excitons and ions
545  G4int NumExcitons =
546  BinomFluct(NumQuanta,fExcitationRatio/(1+fExcitationRatio));
547  G4int NumIons = NumQuanta - NumExcitons;
548 
549  // this section calculates recombination following the modified Birks'
550  // Law of Doke, deposition by deposition, and may be overridden later
551  // in code if a low enough energy necessitates switching to the
552  // Thomas-Imel box model for recombination instead (determined by site)
553  G4double dE, dx=0, LET=0, recombProb;
554  dE = TotalEnergyDeposit/CLHEP::MeV;
555  if ( particleName != "e-" && particleName != "e+" && z1 != z2 &&
556  particleName != "mu-" && particleName != "mu+" ) {
557  //in other words, if it's a gamma,ion,proton,alpha,pion,et al. do not
558  //use the step length provided by Geant4 because it's not relevant,
559  //instead calculate an estimated LET and range of the electrons that
560  //would have been produced if Geant4 could track them
561  LET = CalculateElectronLET( 1000*dE, z1 );
562  if(LET) dx = dE/(Density*LET); //find the range based on the LET
563  if(abs(aParticle->GetPDGcode())==2112) dx=0;
564  }
565  else { //normal case of an e-/+ energy deposition recorded by Geant
566  dx = aStep.GetStepLength()/CLHEP::cm;
567  if(dx) LET = (dE/dx)*(1/Density); //lin. energy xfer (prop. to dE/dx)
568  if ( LET > 0 && dE > 0 && dx > 0 ) {
569  G4double ratio = CalculateElectronLET(dE*1e3,z1)/LET;
570  if ( j == 1 && ratio < 0.7 && !ThomasImelTail &&
571  particleName == "e-" ) {
572  dx /= ratio; LET *= ratio; }}
573  }
574  DokeBirks[1] = DokeBirks[0]/(1-DokeBirks[2]); //B=A/(1-C) (see paper)
575  //Doke/Birks' Law as spelled out in the NEST paper
576  recombProb = (DokeBirks[0]*LET)/(1+DokeBirks[1]*LET)+DokeBirks[2];
577  if ( Phase == kStateLiquid ) {
578  if ( z1 == 54 ) recombProb *= (Density/Density_LXe);
579  if ( z1 == 18 ) recombProb *= (Density/Density_LAr);
580  }
581  //check against unphysicality resulting from rounding errors
582  if(recombProb<0) recombProb=0;
583  if(recombProb>1) recombProb=1;
584  //use binomial distribution to assign photons, electrons, where photons
585  //are excitons plus recombined ionization electrons, while final
586  //collected electrons are the "escape" (non-recombined) electrons
587  G4int NumPhotons = NumExcitons + BinomFluct(NumIons,recombProb);
588  G4int NumElectrons = NumQuanta - NumPhotons;
589 
590  fEnergyDep = TotalEnergyDeposit;
591  fNumIonElectrons = NumElectrons;
592  fNumScintPhotons = NumPhotons;
593 
594  // next section increments the numbers of excitons, ions, photons, and
595  // electrons for the appropriate interaction site; it only appears to
596  // be redundant by saving seemingly no longer needed exciton and ion
597  // counts, these having been already used to calculate the number of ph
598  // and e- above, whereas it does need this later for Thomas-Imel model
599  char numExc[80]; char numIon[80]; char numPho[80]; char numEle[80];
600  sprintf(numExc,"N_EXC_%i",counter); sprintf(numIon,"N_ION_%i",counter);
601  aMaterialPropertiesTable->AddConstProperty( numExc, NumExcitons );
602  aMaterialPropertiesTable->AddConstProperty( numIon, NumIons );
603  sprintf(numPho,"N_PHO_%i",counter); sprintf(numEle,"N_ELE_%i",counter);
604  aMaterialPropertiesTable->AddConstProperty( numPho, NumPhotons );
605  aMaterialPropertiesTable->AddConstProperty( numEle, NumElectrons );
606 
607  // increment and save the total track length, and save interaction
608  // times for later, when generating the scintillation quanta
609  char trackL[80]; char time00[80]; char time01[80]; char energy[80];
610  sprintf(trackL,"TRACK_%i",counter); sprintf(energy,"ENRGY_%i",counter);
611  sprintf(time00,"TIME0_%i",counter); sprintf(time01,"TIME1_%i",counter);
612  delta = aMaterialPropertiesTable->GetConstProperty( trackL );
613  G4double energ = aMaterialPropertiesTable->GetConstProperty( energy );
614  delta += dx*CLHEP::cm; energ += dE*CLHEP::MeV;
615  aMaterialPropertiesTable->AddConstProperty( trackL, delta );
616  aMaterialPropertiesTable->AddConstProperty( energy, energ );
617  if ( TotalEnergyDeposit > 0 ) {
618  G4double deltaTime = aMaterialPropertiesTable->
619  GetConstProperty( time00 );
620  //for charged particles, which continuously lose energy, use initial
621  //interaction time as the minimum time, otherwise use only the final
622  if (aParticle->GetCharge() != 0) {
623  if (t0 < deltaTime)
624  aMaterialPropertiesTable->AddConstProperty( time00, t0 );
625  }
626  else {
627  if (t1 < deltaTime)
628  aMaterialPropertiesTable->AddConstProperty( time00, t1 );
629  }
630  deltaTime = aMaterialPropertiesTable->GetConstProperty( time01 );
631  //find the maximum possible scintillation "birth" time
632  if (t1 > deltaTime)
633  aMaterialPropertiesTable->AddConstProperty( time01, t1 );
634  }
635 
636  // begin the process of setting up creation of scint./ionization
637  TotalEnergyDeposit=aMaterialPropertiesTable->
638  GetConstProperty("ENERGY_DEPOSIT_TOT"); //get the total E deposited
639  InitialKinetEnergy=aMaterialPropertiesTable->
640  GetConstProperty("ENERGY_DEPOSIT_GOL"); //E that should have been
641  if(InitialKinetEnergy > HIENLIM &&
642  abs(aParticle->GetPDGcode()) != 2112) fVeryHighEnergy=true;
643  G4double safety; //margin of error for TotalE.. - InitialKinetEnergy
644  if (fVeryHighEnergy && !fExcitedNucleus) safety = 0.2*CLHEP::keV;
645  else safety = 2.*CLHEP::eV;
646 
647  //force a scintillation dump for NR and for full nuclear de-excitation
648  if( !anExcitationEnergy && pDef->GetParticleType() == "nucleus" &&
649  aTrack.GetTrackStatus() != fAlive && !fAlpha )
650  InitialKinetEnergy = TotalEnergyDeposit;
651  if ( particleName == "neutron" || particleName == "antineutron" )
652  InitialKinetEnergy = TotalEnergyDeposit;
653 
654  //force a dump of all saved scintillation under the following
655  //conditions: energy goal reached, and current particle dead, or an
656  //error has occurred and total has exceeded goal (shouldn't happen)
657  if( std::abs(TotalEnergyDeposit-InitialKinetEnergy)<safety ||
658  TotalEnergyDeposit>=InitialKinetEnergy ){
659  dx = 0; dE = 0;
660  //calculate the total number of quanta from all sites and all
661  //interactions so that the number of secondaries gets set correctly
662  NumPhotons = 0; NumElectrons = 0;
663  for(i=0;i<j;i++) {
664  sprintf(numPho,"N_PHO_%d",i); sprintf(numEle,"N_ELE_%d",i);
665  sprintf(trackL,"TRACK_%d",i); sprintf(energy,"ENRGY_%d",i);
666  //add up track lengths of all sites, for a total LET calc (later)
667  dx += aMaterialPropertiesTable->GetConstProperty(trackL);
668  dE += aMaterialPropertiesTable->GetConstProperty(energy);
669  }
670 
671  G4int buffer = 100; if ( fVeryHighEnergy ) buffer = 1;
672  fParticleChange.SetNumberOfSecondaries(buffer*(NumPhotons+NumElectrons));
673 
674  if (fTrackSecondariesFirst) {
675  if (aTrack.GetTrackStatus() == fAlive )
676  fParticleChange.ProposeTrackStatus(fSuspend);
677  }
678 
679 
680  // begin the loop over all sites which generates all the quanta
681  for(i=0;i<j;i++) {
682  // get the position X,Y,Z, exciton and ion numbers, total track
683  // length of the site, and interaction times
684  sprintf(xCoord,"POS_X_%d",i); sprintf(yCoord,"POS_Y_%d",i);
685  sprintf(zCoord,"POS_Z_%d",i);
686  sprintf(numExc,"N_EXC_%d",i); sprintf(numIon,"N_ION_%d",i);
687  sprintf(numPho,"N_PHO_%d",i); sprintf(numEle,"N_ELE_%d",i);
688  NumExcitons = (G4int)aMaterialPropertiesTable->
689  GetConstProperty( numExc );
690  NumIons = (G4int)aMaterialPropertiesTable->
691  GetConstProperty( numIon );
692  sprintf(trackL,"TRACK_%d",i); sprintf(energy,"ENRGY_%d",i);
693  sprintf(time00,"TIME0_%d",i); sprintf(time01,"TIME1_%d",i);
694  delta = aMaterialPropertiesTable->GetConstProperty( trackL );
695  energ = aMaterialPropertiesTable->GetConstProperty( energy );
696  t0 = aMaterialPropertiesTable->GetConstProperty( time00 );
697  t1 = aMaterialPropertiesTable->GetConstProperty( time01 );
698 
699  //if site is small enough, override the Doke/Birks' model with
700  //Thomas-Imel, but not if we're dealing with super-high energy
701  //particles, and if it's NR force Thomas-Imel (though NR should be
702  //already short enough in track even up to O(100) keV)
703  if ( (delta < R0 && !fVeryHighEnergy) || z2 == z1 || fAlpha ) {
704  if( z1 == 54 && ElectricField && //see NEST paper for ER formula
705  Phase == kStateLiquid ) {
706  if ( abs ( z1 - z2 ) && //electron recoil
707  abs ( aParticle->GetPDGcode() ) != 2112 ) {
708  ThomasImel = 0.056776*pow(ElectricField,-0.11844);
709  if ( fAlpha ) //technically ER, but special
710  ThomasImel=0.057675*pow(ElectricField,-0.49362);
711  } //end electron recoil (ER)
712  else { //nuclear recoil
713  // spline of NR data of C.E. Dahl PhD Thesis Princeton '09
714  // functions found using zunzun.com
715  ThomasImel =
716  -0.15169*pow((ElectricField+215.25),0.01811)+0.20952;
717  } //end NR information
718  // Never let LY exceed 0-field yield!
719  if (ThomasImel > 0.19) ThomasImel = 0.19;
720  if (ThomasImel < 0.00) ThomasImel = 0.00;
721  } //end non-zero E-field segment
722  if ( Phase == kStateLiquid ) {
723  if ( z1 == 54 ) ThomasImel *= pow((Density/2.84),0.3);
724  if ( z1 == 18 ) ThomasImel *= pow((Density/Density_LAr),0.3);
725  }
726  //calculate the Thomas-Imel recombination probability, which
727  //depends on energy via NumIons, but not on dE/dx, and protect
728  //against seg fault by ensuring a positive number of ions
729  if (NumIons > 0) {
730  G4double xi;
731  xi = (G4double(NumIons)/4.)*ThomasImel;
732  if ( InitialKinetEnergy == 9.4*CLHEP::keV ) {
733  G4double NumIonsEff = 1.066e7*pow(t0/CLHEP::ns,-1.303)*
734  (0.17163+162.32/(ElectricField+191.39));
735  if ( NumIonsEff > 1e6 ) NumIonsEff = 1e6;
736  xi = (G4double(NumIonsEff)/4.)*ThomasImel;
737  }
738  if ( fKr83m && ElectricField==0 )
739  xi = (G4double(1.3*NumIons)/4.)*ThomasImel;
740  recombProb = 1-log(1+xi)/xi;
741  if(recombProb<0) recombProb=0;
742  if(recombProb>1) recombProb=1;
743  }
744  //just like Doke: simple binomial distribution
745  NumPhotons = NumExcitons + BinomFluct(NumIons,recombProb);
746  NumElectrons = (NumExcitons + NumIons) - NumPhotons;
747  //override Doke NumPhotons and NumElectrons
748  aMaterialPropertiesTable->
749  AddConstProperty( numPho, NumPhotons );
750  aMaterialPropertiesTable->
751  AddConstProperty( numEle, NumElectrons );
752  }
753 
754  // grab NumPhotons/NumElectrons, which come from Birks if
755  // the Thomas-Imel block of code above was not executed
756  NumPhotons = (G4int)aMaterialPropertiesTable->
757  GetConstProperty( numPho );
758  NumElectrons =(G4int)aMaterialPropertiesTable->
759  GetConstProperty( numEle );
760  // extra Fano factor caused by recomb. fluct.
761  G4double FanoFactor =0; //ionization channel
762  if(Phase == kStateLiquid && fYieldFactor == 1) {
763  FanoFactor =
764  2575.9*pow((ElectricField+15.154),-0.64064)-1.4707;
765  if ( fKr83m ) TotalEnergyDeposit = 4*CLHEP::keV;
766  if ( (dE/CLHEP::keV) <= 100 && ElectricField >= 0 ) {
767  G4double keVee = (TotalEnergyDeposit/(100.*CLHEP::keV));
768  if ( keVee <= 0.06 )
769  FanoFactor *= -0.00075+0.4625*keVee+34.375*pow(keVee,2.);
770  else
771  FanoFactor *= 0.069554+1.7322*keVee-.80215*pow(keVee,2.);
772  }
773  }
774  if ( Phase == kStateGas && Density>0.5 ) FanoFactor =
775  0.42857-4.7857*Density+7.8571*pow(Density,2.);
776  if( FanoFactor <= 0 || fVeryHighEnergy ) FanoFactor = 0;
777  NumQuanta = NumPhotons + NumElectrons;
778  if(z1==54 && FanoFactor) NumElectrons = G4int(
779  floor(GaussGen.fire(NumElectrons,
780  sqrt(FanoFactor*NumElectrons))+0.5));
781  NumPhotons = NumQuanta - NumElectrons;
782  if ( NumElectrons <= 0 ) NumElectrons = 0;
783  if ( NumPhotons <= 0 ) NumPhotons = 0;
784  else { //other effects
785  if ( fAlpha ) //bi-excitonic quenching due to high dE/dx
786  NumPhotons = BinomFluct(NumPhotons,biExc);
787  NumPhotons = BinomFluct(NumPhotons,QE_EFF);
788  } NumElectrons = G4int(floor(NumElectrons*phe_per_e+0.5));
789 
790 
791 
792 
793  // new stuff to make Kr-83m work properly
794  if(fKr83m || InitialKinetEnergy==9.4*CLHEP::keV) fKr83m += dE/CLHEP::keV;
795  if(fKr83m > 41) fKr83m = 0;
796  if ( SinglePhase ) //for a 1-phase det. don't propagate e-'s
797  NumElectrons = 0; //saves simulation time
798 
799  // reset material properties numExc, numIon, numPho, numEle, as
800  // their values have been used or stored elsewhere already
801  aMaterialPropertiesTable->AddConstProperty( numExc, 0 );
802  aMaterialPropertiesTable->AddConstProperty( numIon, 0 );
803  aMaterialPropertiesTable->AddConstProperty( numPho, 0 );
804  aMaterialPropertiesTable->AddConstProperty( numEle, 0 );
805 
806  // start particle creation loop
807  if( InitialKinetEnergy < MAX_ENE && InitialKinetEnergy > MIN_ENE &&
808  !fMultipleScattering )
809  NumQuanta = NumPhotons + NumElectrons;
810  else NumQuanta = 0;
811  for(k = 0; k < NumQuanta; k++) {
812  G4double sampledEnergy;
813  std::unique_ptr<G4DynamicParticle> aQuantum;
814 
815  // Generate random direction
816  G4double cost = 1. - 2.*UniformGen.fire();
817  G4double sint = std::sqrt((1.-cost)*(1.+cost));
818  G4double phi = CLHEP::twopi*UniformGen.fire();
819  G4double sinp = std::sin(phi); G4double cosp = std::cos(phi);
820  G4double px = sint*cosp; G4double py = sint*sinp;
821  G4double pz = cost;
822 
823  // Create momentum direction vector
824  G4ParticleMomentum photonMomentum(px, py, pz);
825 
826  // case of photon-specific stuff
827  if (k < NumPhotons) {
828  // Determine polarization of new photon
829  G4double sx = cost*cosp;
830  G4double sy = cost*sinp;
831  G4double sz = -sint;
832  G4ThreeVector photonPolarization(sx, sy, sz);
833  G4ThreeVector perp = photonMomentum.cross(photonPolarization);
834  phi = CLHEP::twopi*UniformGen.fire();
835  sinp = std::sin(phi);
836  cosp = std::cos(phi);
837  photonPolarization = cosp * photonPolarization + sinp * perp;
838  photonPolarization = photonPolarization.unit();
839 
840  // Generate a new photon or electron:
841  sampledEnergy = GaussGen.fire(PhotMean,PhotWidth);
842  aQuantum = std::make_unique<G4DynamicParticle>(G4OpticalPhoton::OpticalPhoton(),
843  photonMomentum);
844  aQuantum->SetPolarization(photonPolarization.x(),
845  photonPolarization.y(),
846  photonPolarization.z());
847  }
848 
849  else { // this else statement is for ionization electrons
850  if(ElectricField) {
851  // point all electrons straight up, for drifting
852  G4ParticleMomentum electronMomentum(0, 0, -FieldSign);
853  aQuantum = std::make_unique<G4DynamicParticle>(G4ThermalElectron::ThermalElectron(),
854  electronMomentum);
855  if ( Phase == kStateGas ) {
856  sampledEnergy = GetGasElectronDriftSpeed(ElectricField,nDensity);
857  }
858  else
859  sampledEnergy = GetLiquidElectronDriftSpeed(Temperature,ElectricField,MillerDriftSpeed,z1);
860  }
861  else {
862  // use "photonMomentum" for the electrons in the case of zero
863  // electric field, which is just randomized vector we made
864  aQuantum = std::make_unique<G4DynamicParticle>(G4ThermalElectron::ThermalElectron(),
865  photonMomentum);
866  sampledEnergy = 1.38e-23*(CLHEP::joule/CLHEP::kelvin)*Temperature;
867  }
868  }
869 
870  //assign energy to make particle real
871  aQuantum->SetKineticEnergy(sampledEnergy);
872 
873  // Generate new G4Track object:
874  // emission time distribution
875 
876  // first an initial birth time is provided that is typically
877  // <<1 ns after the initial interaction in the simulation, then
878  // singlet, triplet lifetimes, and recombination time, are
879  // handled here, to create a realistic S1 pulse shape/timing
880  G4double aSecondaryTime = t0+UniformGen.fire()*(t1-t0)+evtStrt;
881  if (tau1<0) { tau1=0; }
882  if (tau3<0) { tau3=0; }
883  if (tauR<0) { tauR=0; }
884  if ( aQuantum->GetDefinition()->
885  GetParticleName()=="opticalphoton" ) {
886  if ( abs(z2-z1) && !fAlpha && //electron recoil
887  abs(aParticle->GetPDGcode()) != 2112 ) {
888  LET = (energ/CLHEP::MeV)/(delta/CLHEP::cm)*(1/Density); //avg LET over all
889  //in future, this will be done interaction by interaction
890  // Next, find the recombination time, which is LET-dependent
891  // via ionization density (Kubota et al. Phys. Rev. B 20
892  // (1979) 3486). We find the LET-dependence by fitting to the
893  // E-dependence (Akimov et al. Phys. Lett. B 524 (2002) 245).
894  if ( Phase == kStateLiquid && z1 == 54 )
895  tauR = 3.5*((1+0.41*LET)/(0.18*LET))*CLHEP::ns
896  *exp(-0.00900*ElectricField);
897  //field dependence based on fitting Fig. 9 of Dawson et al.
898  //NIM A 545 (2005) 690
899  //singlet-triplet ratios adapted from Kubota 1979, converted
900  //into correct units to work here, and separately done for
901  //excitation and recombination processes for electron recoils
902  //and assumed same for all LET (may vary)
903  SingTripRatioX = GaussGen.fire(0.17,0.05);
904  SingTripRatioR = GaussGen.fire(0.8,0.2);
905  if ( z1 == 18 ) {
906  SingTripRatioR = 0.2701+0.003379*LET-4.7338e-5*pow(LET,2.)
907  +8.1449e-6*pow(LET,3.); SingTripRatioX = SingTripRatioR;
908  if( LET < 3 ) {
909  SingTripRatioX = GaussGen.fire(0.36,0.06);
910  SingTripRatioR = GaussGen.fire(0.5,0.2); }
911  }
912  }
913  else if ( fAlpha ) { //alpha particles
914  SingTripRatioR = GaussGen.fire(2.3,0.51);
915  //currently based on Dawson 05 and Tey. 11 (arXiv:1103.3689)
916  //real ratio is likely a gentle function of LET
917  if (z1==18) SingTripRatioR = (-0.065492+1.9996
918  *exp(-dE/CLHEP::MeV))/(1+0.082154/pow(dE/CLHEP::MeV,2.)) + 2.1811;
919  SingTripRatioX = SingTripRatioR;
920  }
921  else { //nuclear recoil
922  //based loosely on Hitachi et al. Phys. Rev. B 27 (1983) 5279
923  //with an eye to reproducing Akimov 2002 Fig. 9
924  SingTripRatioR = GaussGen.fire(7.8,1.5);
925  if (z1==18) SingTripRatioR = 0.22218*pow(energ/CLHEP::keV,0.48211);
926  SingTripRatioX = SingTripRatioR;
927  }
928  // now, use binomial distributions to determine singlet and
929  // triplet states (and do separately for initially excited guys
930  // and recombining)
931  if ( k > NumExcitons ) {
932  //the recombination time is non-exponential, but approximates
933  //to exp at long timescales (see Kubota '79)
934  aSecondaryTime += tauR*(1./UniformGen.fire()-1);
935  if(UniformGen.fire()<SingTripRatioR/(1+SingTripRatioR))
936  aSecondaryTime -= tau1*log(UniformGen.fire());
937  else aSecondaryTime -= tau3*log(UniformGen.fire());
938  }
939  else {
940  if(UniformGen.fire()<SingTripRatioX/(1+SingTripRatioX))
941  aSecondaryTime -= tau1*log(UniformGen.fire());
942  else aSecondaryTime -= tau3*log(UniformGen.fire());
943  }
944  }
945  else { //electron trapping at the liquid/gas interface
946  G4double gainField = 12;
947  G4double tauTrap = 884.83-62.069*gainField;
948  if ( Phase == kStateLiquid )
949  aSecondaryTime -= tauTrap*CLHEP::ns*log(UniformGen.fire());
950  }
951 
952  // emission position distribution --
953  // Generate the position of a new photon or electron, with NO
954  // stochastic variation because that could lead to particles
955  // being mistakenly generated outside of your active region by
956  // Geant4, but real-life finite detector position resolution
957  // wipes out any effects from here anyway...
958  std::string sx(xCoord);
959  std::string sy(yCoord);
960  std::string sz(zCoord);
961  x0[0] = aMaterialPropertiesTable->GetConstProperty( xCoord );
962  x0[1] = aMaterialPropertiesTable->GetConstProperty( yCoord );
963  x0[2] = aMaterialPropertiesTable->GetConstProperty( zCoord );
964  G4double radius = sqrt(pow(x0[0],2.)+pow(x0[1],2.));
965  //re-scale radius to ensure no generation of quanta outside
966  //the active volume of your simulation due to Geant4 rounding
967  if ( radius >= R_TOL ) {
968  if (x0[0] == 0) { x0[0] = 1*CLHEP::nm; }
969  if (x0[1] == 0) { x0[1] = 1*CLHEP::nm; }
970  radius -= R_TOL; phi = atan ( x0[1] / x0[0] );
971  x0[0] = fabs(radius*cos(phi))*((fabs(x0[0]))/(x0[0]));
972  x0[1] = fabs(radius*sin(phi))*((fabs(x0[1]))/(x0[1]));
973  }
974  //position of the new secondary particle is ready for use
975  G4ThreeVector aSecondaryPosition = x0;
976  if ( k >= NumPhotons && diffusion && ElectricField > 0 ) {
977  G4double D_T = 64*pow(1e-3*ElectricField,-.17);
978  //fit to Aprile and Doke 2009, arXiv:0910.4956 (Fig. 12)
979  G4double D_L = 13.859*pow(1e-3*ElectricField,-0.58559);
980  //fit to Aprile and Doke and Sorensen 2011, arXiv:1102.2865
981  if ( Phase == kStateLiquid && z1 == 18 ) {
982  D_T = 93.342*pow(ElectricField/nDensity,0.041322);
983  D_L = 0.15 * D_T; }
984  if ( Phase == kStateGas && z1 == 54 ) {
985  D_L=4.265+19097/ElectricField-1.7397e6/pow(ElectricField,2.)+
986  1.2477e8/pow(ElectricField,3.); D_T *= 0.01;
987  } D_T *= CLHEP::cm2/CLHEP::s; D_L *= CLHEP::cm2/CLHEP::s;
988  if (ElectricField<100 && Phase == kStateLiquid) D_L = 8*CLHEP::cm2/CLHEP::s;
989  G4double vDrift = sqrt((2*sampledEnergy)/(EMASS));
990  if ( BORDER == 0 ) x0[2] = 0;
991  G4double sigmaDT = sqrt(2*D_T*fabs(BORDER-x0[2])/vDrift);
992  G4double sigmaDL = sqrt(2*D_L*fabs(BORDER-x0[2])/vDrift);
993  G4double dr = std::abs(GaussGen.fire(0.,sigmaDT));
994  phi = CLHEP::twopi * UniformGen.fire();
995  aSecondaryPosition[0] += cos(phi) * dr;
996  aSecondaryPosition[1] += sin(phi) * dr;
997  aSecondaryPosition[2] += GaussGen.fire(0.,sigmaDL);
998  radius = std::sqrt(std::pow(aSecondaryPosition[0],2.)+
999  std::pow(aSecondaryPosition[1],2.));
1000  if(aSecondaryPosition[2] >= BORDER && Phase == kStateLiquid) {
1001  if ( BORDER != 0 ) aSecondaryPosition[2] = BORDER - R_TOL; }
1002  if(aSecondaryPosition[2] <= PMT && !GlobalFields)
1003  aSecondaryPosition[2] = PMT + R_TOL;
1004  } //end of electron diffusion code
1005 
1006  // GEANT4 business: stuff you need to make a new track
1007  if ( aSecondaryTime < 0 ) aSecondaryTime = 0; //no neg. time
1008  /*
1009  auto aSecondaryTrack =
1010  std::make_unique<G4Track>(aQuantum,aSecondaryTime,aSecondaryPosition);
1011  if ( k < NumPhotons || radius < R_MAX )
1012  {
1013 
1014  */
1015  }
1016 
1017  //reset bunch of things when done with an interaction site
1018  aMaterialPropertiesTable->AddConstProperty( xCoord, 999*CLHEP::km );
1019  aMaterialPropertiesTable->AddConstProperty( yCoord, 999*CLHEP::km );
1020  aMaterialPropertiesTable->AddConstProperty( zCoord, 999*CLHEP::km );
1021  aMaterialPropertiesTable->AddConstProperty( trackL, 0*CLHEP::um );
1022  aMaterialPropertiesTable->AddConstProperty( energy, 0*CLHEP::eV );
1023  aMaterialPropertiesTable->AddConstProperty( time00, DBL_MAX );
1024  aMaterialPropertiesTable->AddConstProperty( time01, -1*CLHEP::ns );
1025 
1026  } //end of interaction site loop
1027 
1028  //more things to reset...
1029  aMaterialPropertiesTable->
1030  AddConstProperty( "TOTALNUM_INT_SITES", 0 );
1031  aMaterialPropertiesTable->
1032  AddConstProperty( "ENERGY_DEPOSIT_TOT", 0*CLHEP::keV );
1033  aMaterialPropertiesTable->
1034  AddConstProperty( "ENERGY_DEPOSIT_GOL", 0*CLHEP::MeV );
1035  fExcitedNucleus = false;
1036  fAlpha = false;
1037  }
1038 
1039  //don't do anything when you're not ready to scintillate
1040  else {
1041  fParticleChange.SetNumberOfSecondaries(0);
1042  return fParticleChange;
1043  }
1044 
1045  return fParticleChange;
1046 }
t0
code to link reconstructed objects back to the MC truth information
Definition: DirectHitParticleAssns_tool.cc:18
s
Float_t s
Definition: plot.C:23
AVO
#define AVO
Definition: NestAlg.cxx:1
Density_LAr
#define Density_LAr
Definition: NestAlg.cxx:26
ThomasImelTail
G4bool ThomasImelTail
Definition: NestAlg.cxx:43
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:1060
x1
Float_t x1[n_points_granero]
Definition: compare.C:5
biExc
G4double biExc
Definition: NestAlg.cxx:45
SinglePhase
G4bool SinglePhase
Definition: NestAlg.cxx:43
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:37
NestAlg::fParticleChange
G4VParticleChange fParticleChange
pointer to G4VParticleChange
Definition: NestAlg.h:43
NestAlg::fNumIonElectrons
int fNumIonElectrons
number of ionization electrons produced by step
Definition: NestAlg.h:41
NestAlg::BinomFluct
G4int BinomFluct(G4int N0, G4double prob)
Definition: NestAlg.cxx:1190
PMT
#define PMT
Definition: NestAlg.cxx:19
NestAlg::fNumScintPhotons
int fNumScintPhotons
number of photons produced by the step
Definition: NestAlg.h:40
a2
#define a2
Definition: HoughBaseAlg.cxx:67
G4ThermalElectron::ThermalElectron
static G4ThermalElectron * ThermalElectron()
Definition: G4ThermalElectron.cxx:92
diffusion
G4bool diffusion
Definition: NestAlg.cxx:41
energy
double energy
Definition: plottest35.C:25
Density_LXe
#define Density_LXe
Definition: NestAlg.cxx:25
radius
Double_t radius
Definition: TG43_relative_dose.C:19
NestAlg::GetGasElectronDriftSpeed
G4double GetGasElectronDriftSpeed(G4double efieldinput, G4double density)
Definition: NestAlg.cxx:1049
NestAlg::fEngine
CLHEP::HepRandomEngine & fEngine
random engine
Definition: NestAlg.h:47
NestAlg::CalculateElectronLET
G4double CalculateElectronLET(G4double E, G4int Z)
Definition: NestAlg.cxx:1162
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:1214
BOT
#define BOT
Definition: NestAlg.cxx:18
MillerDriftSpeed
#define MillerDriftSpeed
Definition: NestAlg.cxx:3
NestAlg::fElementPropInit
std::map< int, bool > fElementPropInit
Definition: NestAlg.h:44
ANE
#define ANE
Definition: NestAlg.cxx:14
NestAlg::fExcitationRatio
double fExcitationRatio
Definition: NestAlg.h:38
CTH
#define CTH
Definition: NestAlg.cxx:17
e
Float_t e
Definition: plot.C:34
OutElectrons
G4bool OutElectrons
Definition: NestAlg.cxx:43
TOP
#define TOP
Definition: NestAlg.cxx:13
MIN_ENE
#define MIN_ENE
Definition: NestAlg.cxx:20
NestAlg::fEnergyDep
double fEnergyDep
energy deposited by the step
Definition: NestAlg.h:42
a1
#define a1
Definition: HoughBaseAlg.cxx:66
double NestAlg::EnergyDeposition ( ) const
inline

Definition at line 26 of file NestAlg.h.

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

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

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

Definition at line 1049 of file NestAlg.cxx.

Referenced by CalculateIonizationAndScintillation().

1051 {
1052  std::cout << "WARNING: NestAlg::GetGasElectronDriftSpeed(G4double, G4double) "
1053  << "is not defined, returning bogus value of -999." << std::endl;
1054 
1055  return -999.;
1056 }
G4double NestAlg::GetLiquidElectronDriftSpeed ( double  T,
double  F,
G4bool  M,
G4int  Z 
)
private

Definition at line 1060 of file NestAlg.cxx.

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

Referenced by CalculateIonizationAndScintillation().

1063  {
1064  if(efieldinput<0) efieldinput *= (-1);
1065  //Liquid equation one (165K) coefficients
1066  G4double onea=144623.235704015,
1067  oneb=850.812714257629,
1068  onec=1192.87056676815,
1069  oned=-395969.575204061,
1070  onef=-355.484170008875,
1071  oneg=-227.266219627672,
1072  oneh=223831.601257495,
1073  onei=6.1778950907965,
1074  onej=18.7831533426398,
1075  onek=-76132.6018884368;
1076  //Liquid equation two (200K) coefficients
1077  G4double twoa=17486639.7118995,
1078  twob=-113.174284723134,
1079  twoc=28.005913193763,
1080  twod=167994210.094027,
1081  twof=-6766.42962575088,
1082  twog=901.474643115395,
1083  twoh=-185240292.471665,
1084  twoi=-633.297790813084,
1085  twoj=87.1756135457949;
1086  //Liquid equation three (230K) coefficients
1087  G4double thra=10626463726.9833,
1088  thrb=224025158.134792,
1089  thrc=123254826.300172,
1090  thrd=-4563.5678061122,
1091  thrf=-1715.269592063,
1092  thrg=-694181.921834368,
1093  thrh=-50.9753281079838,
1094  thri=58.3785811395493,
1095  thrj=201512.080026704;
1096  G4double y1=0,y2=0,f1=0,f2=0,f3=0,edrift=0,
1097  t1=0,t2=0,slope=0,intercept=0;
1098 
1099  //Equations defined
1100  f1=onea/(1+exp(-(efieldinput-oneb)/onec))+oned/
1101  (1+exp(-(efieldinput-onef)/oneg))+
1102  oneh/(1+exp(-(efieldinput-onei)/onej))+onek;
1103  f2=twoa/(1+exp(-(efieldinput-twob)/twoc))+twod/
1104  (1+exp(-(efieldinput-twof)/twog))+
1105  twoh/(1+exp(-(efieldinput-twoi)/twoj));
1106  f3=thra*exp(-thrb*efieldinput)+thrc*exp(-(pow(efieldinput-thrd,2))/
1107  (thrf*thrf))+
1108  thrg*exp(-(pow(efieldinput-thrh,2)/(thri*thri)))+thrj;
1109 
1110  if(efieldinput<20 && efieldinput>=0) {
1111  f1=2951*efieldinput;
1112  f2=5312*efieldinput;
1113  f3=7101*efieldinput;
1114  }
1115  //Cases for tempinput decides which 2 equations to use lin. interpolation
1116  if(tempinput<200.0 && tempinput>165.0) {
1117  y1=f1;
1118  y2=f2;
1119  t1=165.0;
1120  t2=200.0;
1121  }
1122  if(tempinput<230.0 && tempinput>200.0) {
1123  y1=f2;
1124  y2=f3;
1125  t1=200.0;
1126  t2=230.0;
1127  }
1128  if((tempinput>230.0 || tempinput<165.0) && !Miller) {
1129  G4cout << "\nWARNING: TEMPERATURE OUT OF RANGE (165-230 K)\n";
1130  return 0;
1131  }
1132  if (tempinput == 165.0) edrift = f1;
1133  else if (tempinput == 200.0) edrift = f2;
1134  else if (tempinput == 230.0) edrift = f3;
1135  else { //Linear interpolation
1136  //frac=(tempinput-t1)/(t2-t1);
1137  slope = (y1-y2)/(t1-t2);
1138  intercept=y1-slope*t1;
1139  edrift=slope*tempinput+intercept;
1140  }
1141 
1142  if ( Miller ) {
1143  if ( efieldinput <= 40. )
1144  edrift = -0.13274+0.041082*efieldinput-0.0006886*pow(efieldinput,2.)+
1145  5.5503e-6*pow(efieldinput,3.);
1146  else
1147  edrift = 0.060774*efieldinput/pow(1+0.11336*pow(efieldinput,0.5218),2.);
1148  if ( efieldinput >= 1e5 ) edrift = 2.7;
1149  if ( efieldinput >= 100 )
1150  edrift -= 0.017 * ( tempinput - 163 );
1151  else
1152  edrift += 0.017 * ( tempinput - 163 );
1153  edrift *= 1e5; //put into units of cm/sec. from mm/usec.
1154  }
1155  if ( Z == 18 ) edrift = 1e5 * (.097384*pow(log10(efieldinput),3.0622)-.018614*sqrt(efieldinput) );
1156  if ( edrift < 0 ) edrift = 0.;
1157  edrift = 0.5*EMASS*pow(edrift*CLHEP::cm/CLHEP::s,2.);
1158  return edrift;
1159 }
s
Float_t s
Definition: plot.C:23
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
larg4::f3
DecayPhysicsFactory f3
Definition: CustomPhysicsBuiltIns.cxx:12
larg4::f2
SynchrotronAndGN f2
Definition: CustomPhysicsBuiltIns.cxx:11
y2
Float_t y2[n_points_geant4]
Definition: compare.C:26
Z
Float_t Z
Definition: plot.C:39
larg4::f1
EmPhysicsFactory f1
Definition: CustomPhysicsBuiltIns.cxx:10
t2
TTree * t2
Definition: plottest35.C:36
void NestAlg::InitMatPropValues ( G4MaterialPropertiesTable *  nobleElementMat,
int  z 
)
private

Definition at line 1214 of file NestAlg.cxx.

References energy, fElementPropInit, and z.

Referenced by CalculateIonizationAndScintillation().

1215  {
1216  char xCoord[80]; char yCoord[80]; char zCoord[80];
1217  char numExc[80]; char numIon[80]; char numPho[80]; char numEle[80];
1218  char trackL[80]; char time00[80]; char time01[80]; char energy[80];
1219 
1220  // for loop to initialize the interaction site mat'l properties
1221  for( G4int i=0; i<10000; i++ ) {
1222  sprintf(xCoord,"POS_X_%d",i); sprintf(yCoord,"POS_Y_%d",i);
1223  sprintf(zCoord,"POS_Z_%d",i);
1224  nobleElementMat->AddConstProperty( xCoord, 999*CLHEP::km );
1225  nobleElementMat->AddConstProperty( yCoord, 999*CLHEP::km );
1226  nobleElementMat->AddConstProperty( zCoord, 999*CLHEP::km );
1227  sprintf(numExc,"N_EXC_%d",i); sprintf(numIon,"N_ION_%d",i);
1228  sprintf(numPho,"N_PHO_%d",i); sprintf(numEle,"N_ELE_%d",i);
1229  nobleElementMat->AddConstProperty( numExc, 0 );
1230  nobleElementMat->AddConstProperty( numIon, 0 );
1231  nobleElementMat->AddConstProperty( numPho, 0 );
1232  nobleElementMat->AddConstProperty( numEle, 0 );
1233  sprintf(trackL,"TRACK_%d",i); sprintf(energy,"ENRGY_%d",i);
1234  sprintf(time00,"TIME0_%d",i); sprintf(time01,"TIME1_%d",i);
1235  nobleElementMat->AddConstProperty( trackL, 0*CLHEP::um );
1236  nobleElementMat->AddConstProperty( energy, 0*CLHEP::eV );
1237  nobleElementMat->AddConstProperty( time00, DBL_MAX );
1238  nobleElementMat->AddConstProperty( time01,-1*CLHEP::ns );
1239  }
1240 
1241  // we initialize the total number of interaction sites, a variable for
1242  // updating the amount of energy deposited thus far in the medium, and a
1243  // variable for storing the amount of energy expected to be deposited
1244  nobleElementMat->AddConstProperty( "TOTALNUM_INT_SITES", 0 );
1245  nobleElementMat->AddConstProperty( "ENERGY_DEPOSIT_TOT", 0*CLHEP::keV );
1246  nobleElementMat->AddConstProperty( "ENERGY_DEPOSIT_GOL", 0*CLHEP::MeV );
1247 
1248  fElementPropInit[z] = true;
1249 
1250  return;
1251 }
z
Double_t z
Definition: plot.C:279
energy
double energy
Definition: plottest35.C:25
NestAlg::fElementPropInit
std::map< int, bool > fElementPropInit
Definition: NestAlg.h:44
int NestAlg::NumberIonizationElectrons ( ) const
inline

Definition at line 25 of file NestAlg.h.

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

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

Definition at line 24 of file NestAlg.h.

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

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

Definition at line 23 of file NestAlg.h.

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

Definition at line 22 of file NestAlg.h.

22 { fYieldFactor = yf; }
NestAlg::fYieldFactor
double fYieldFactor
turns scint. on/off
Definition: NestAlg.h:37
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 44 of file NestAlg.h.

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

double NestAlg::fEnergyDep
private

energy deposited by the step

Definition at line 42 of file NestAlg.h.

Referenced by CalculateIonizationAndScintillation().

CLHEP::HepRandomEngine& NestAlg::fEngine
private

random engine

Definition at line 47 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 38 of file NestAlg.h.

Referenced by CalculateIonizationAndScintillation().

int NestAlg::fNumIonElectrons
private

number of ionization electrons produced by step

Definition at line 41 of file NestAlg.h.

Referenced by CalculateIonizationAndScintillation().

int NestAlg::fNumScintPhotons
private

number of photons produced by the step

Definition at line 40 of file NestAlg.h.

Referenced by CalculateIonizationAndScintillation().

G4VParticleChange NestAlg::fParticleChange
private

pointer to G4VParticleChange

Definition at line 43 of file NestAlg.h.

Referenced by CalculateIonizationAndScintillation().

double NestAlg::fYieldFactor
private

turns scint. on/off

Definition at line 37 of file NestAlg.h.

Referenced by CalculateIonizationAndScintillation().

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