Public Member Functions | Private Member Functions | Private Attributes | List of all members
gar::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
 excitons to ions More...
 
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 13 of file NestAlg.h.

Constructor & Destructor Documentation

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

Definition at line 52 of file NestAlg.cxx.

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

Definition at line 68 of file NestAlg.cxx.

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

Member Function Documentation

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

Definition at line 1194 of file NestAlg.cxx.

1194  {
1195  CLHEP::RandGauss GaussGen(fEngine);
1196  CLHEP::RandFlat UniformGen(fEngine);
1197 
1198  G4double mean = N0*prob;
1199  G4double sigma = sqrt(N0*prob*(1-prob));
1200  G4int N1 = 0;
1201  if ( prob == 0.00 ) return N1;
1202  if ( prob == 1.00 ) return N0;
1203 
1204  if ( N0 < 10 ) {
1205  for(G4int i = 0; i < N0; i++) {
1206  if(UniformGen.fire() < prob) N1++;
1207  }
1208  }
1209  else {
1210  N1 = G4int(floor(GaussGen.fire(mean,sigma)+0.5));
1211  }
1212  if ( N1 > N0 ) N1 = N0;
1213  if ( N1 < 0 ) N1 = 0;
1214  return N1;
1215  }
gar::NestAlg::fEngine
CLHEP::HepRandomEngine & fEngine
random engine
Definition: NestAlg.h:49
N1
Definition: group.cpp:53
cet::sqlite::mean
double mean(sqlite3 *db, std::string const &table_name, std::string const &column_name)
Definition: statistics.cc:16
G4double NestAlg::CalculateElectronLET ( G4double  E,
G4int  Z 
)
private

Definition at line 1166 of file NestAlg.cxx.

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

Definition at line 84 of file NestAlg.cxx.

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

Definition at line 28 of file NestAlg.h.

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

Definition at line 1053 of file NestAlg.cxx.

1055  {
1056  std::cout << "WARNING: NestAlg::GetGasElectronDriftSpeed(G4double, G4double) "
1057  << "is not defined, returning bogus value of -999." << std::endl;
1058 
1059  return -999.;
1060  }
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
G4double gar::NestAlg::GetLiquidElectronDriftSpeed ( double  T,
double  F,
G4bool  M,
G4int  Z 
)
private
void NestAlg::InitMatPropValues ( G4MaterialPropertiesTable *  nobleElementMat,
int  z 
)
private

Definition at line 1218 of file NestAlg.cxx.

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

Definition at line 27 of file NestAlg.h.

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

Definition at line 26 of file NestAlg.h.

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

Definition at line 25 of file NestAlg.h.

25 { fExcitationRatio = er; }
gar::NestAlg::fExcitationRatio
double fExcitationRatio
excitons to ions
Definition: NestAlg.h:40
void gar::NestAlg::SetScintillationYieldFactor ( double const &  yf)
inline

Definition at line 24 of file NestAlg.h.

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

Member Data Documentation

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

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

Definition at line 46 of file NestAlg.h.

double gar::NestAlg::fEnergyDep
private

energy deposited by the step

Definition at line 44 of file NestAlg.h.

CLHEP::HepRandomEngine& gar::NestAlg::fEngine
private

random engine

Definition at line 49 of file NestAlg.h.

double gar::NestAlg::fExcitationRatio
private

excitons to ions

N_ex/N_i, the dimensionless ratio of initial

Definition at line 40 of file NestAlg.h.

int gar::NestAlg::fNumIonElectrons
private

number of ionization electrons produced by step

Definition at line 43 of file NestAlg.h.

int gar::NestAlg::fNumScintPhotons
private

number of photons produced by the step

Definition at line 42 of file NestAlg.h.

G4VParticleChange gar::NestAlg::fParticleChange
private

pointer to G4VParticleChange

Definition at line 45 of file NestAlg.h.

double gar::NestAlg::fYieldFactor
private

turns scint. on/off

Definition at line 39 of file NestAlg.h.

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.11