Public Member Functions | Protected Attributes | List of all members
G4S1Light Class Reference

#include <G4S1Light.hh>

Inheritance diagram for G4S1Light:

Public Member Functions

 G4S1Light (const G4String &processName="S1", G4ProcessType type=fElectromagnetic)
 
 ~G4S1Light ()
 
G4bool IsApplicable (const G4ParticleDefinition &aParticleType)
 
G4double GetMeanFreePath (const G4Track &aTrack, G4double, G4ForceCondition *)
 
G4double GetMeanLifeTime (const G4Track &aTrack, G4ForceCondition *)
 
G4VParticleChange * PostStepDoIt (const G4Track &aTrack, const G4Step &aStep)
 
G4VParticleChange * AtRestDoIt (const G4Track &aTrack, const G4Step &aStep)
 
void SetTrackSecondariesFirst (const G4bool state)
 
G4bool GetTrackSecondariesFirst () const
 
void SetScintillationYieldFactor (const G4double yieldfactor)
 
G4double GetScintillationYieldFactor () const
 
void SetScintillationExcitationRatio (const G4double excitationratio)
 
G4double GetScintillationExcitationRatio () const
 

Protected Attributes

G4bool fTrackSecondariesFirst
 
G4bool fExcitedNucleus
 
G4bool fAlpha
 
G4bool fVeryHighEnergy
 
G4bool fMultipleScattering
 
G4double fKr83m
 
G4double YieldFactor
 
G4double ExcitationRatio
 

Detailed Description

Definition at line 45 of file G4S1Light.hh.

Constructor & Destructor Documentation

G4S1Light::G4S1Light ( const G4String &  processName = "S1",
G4ProcessType  type = fElectromagnetic 
)

Definition at line 82 of file G4S1Light.cc.

84  : G4VRestDiscreteProcess(processName, type)
85 {
86  //luxManager = LUXSimManager::GetManager();
87  SetProcessSubType(fScintillation);
88 
89  fTrackSecondariesFirst = false;
90  //particles die first, then scintillation is generated
91 
92  if (verboseLevel>0) {
93  G4cout << GetProcessName() << " is created " << G4endl;
94  }
95 }
G4S1Light::fTrackSecondariesFirst
G4bool fTrackSecondariesFirst
Definition: G4S1Light.hh:109
G4S1Light::~G4S1Light ( )

Definition at line 97 of file G4S1Light.cc.

97 {} //destructor needed to avoid linker error

Member Function Documentation

G4VParticleChange * G4S1Light::AtRestDoIt ( const G4Track &  aTrack,
const G4Step &  aStep 
)

Definition at line 100 of file G4S1Light.cc.

104 {
105  return G4S1Light::PostStepDoIt(aTrack, aStep);
106 }
G4S1Light::PostStepDoIt
G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep)
Definition: G4S1Light.cc:109
G4double G4S1Light::GetMeanFreePath ( const G4Track &  aTrack,
G4double  ,
G4ForceCondition *  condition 
)

Definition at line 1094 of file G4S1Light.cc.

1097 {
1098  *condition = StronglyForced;
1099  // what this does is enforce the G4S1Light physics process as always
1100  // happening, so in effect scintillation is a meta-process on top of
1101  // any and all other energy depositions which may occur, just like the
1102  // original G4Scintillation (disregard DBL_MAX, this function makes the
1103  // mean free path zero really, not infinite)
1104 
1105  return DBL_MAX; //a C-defined constant
1106 }
G4double G4S1Light::GetMeanLifeTime ( const G4Track &  aTrack,
G4ForceCondition *  condition 
)

Definition at line 1110 of file G4S1Light.cc.

1112 {
1113  *condition = Forced;
1114  // this function and this condition has the same effect as the above
1115  return DBL_MAX;
1116 }
G4double G4S1Light::GetScintillationExcitationRatio ( ) const
inline

Definition at line 168 of file G4S1Light.hh.

169 {
170  return ExcitationRatio;
171 }
G4S1Light::ExcitationRatio
G4double ExcitationRatio
Definition: G4S1Light.hh:114
G4double G4S1Light::GetScintillationYieldFactor ( ) const
inline

Definition at line 156 of file G4S1Light.hh.

157 {
158  return YieldFactor;
159 }
G4S1Light::YieldFactor
G4double YieldFactor
Definition: G4S1Light.hh:113
G4bool G4S1Light::GetTrackSecondariesFirst ( ) const
inline

Definition at line 144 of file G4S1Light.hh.

145 {
146  return fTrackSecondariesFirst;
147 }
G4S1Light::fTrackSecondariesFirst
G4bool fTrackSecondariesFirst
Definition: G4S1Light.hh:109
G4bool G4S1Light::IsApplicable ( const G4ParticleDefinition &  aParticleType)
inline

Definition at line 124 of file G4S1Light.hh.

125 {
126  if (aParticleType.GetParticleName() == "opticalphoton") return false;
127  if (aParticleType.IsShortLived()) return false;
128  if (aParticleType.GetParticleName() == "thermalelectron") return false;
129  //if(abs(aParticleType.GetPDGEncoding())==2112 || //neutron (no E-dep.)
130  if(abs(aParticleType.GetPDGEncoding())==12 || //neutrinos (ditto)
131  abs(aParticleType.GetPDGEncoding())==14 ||
132  abs(aParticleType.GetPDGEncoding())==16) return false;
133 
134  return true;
135 }
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
G4VParticleChange * G4S1Light::PostStepDoIt ( const G4Track &  aTrack,
const G4Step &  aStep 
)

Definition at line 109 of file G4S1Light.cc.

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

Definition at line 162 of file G4S1Light.hh.

163 {
164  ExcitationRatio = excitationratio;
165 }
G4S1Light::ExcitationRatio
G4double ExcitationRatio
Definition: G4S1Light.hh:114
void G4S1Light::SetScintillationYieldFactor ( const G4double  yieldfactor)
inline

Definition at line 150 of file G4S1Light.hh.

151 {
152  YieldFactor = yieldfactor;
153 }
G4S1Light::YieldFactor
G4double YieldFactor
Definition: G4S1Light.hh:113
void G4S1Light::SetTrackSecondariesFirst ( const G4bool  state)
inline

Definition at line 138 of file G4S1Light.hh.

139 {
140  fTrackSecondariesFirst = state;
141 }
G4S1Light::fTrackSecondariesFirst
G4bool fTrackSecondariesFirst
Definition: G4S1Light.hh:109

Member Data Documentation

G4double G4S1Light::ExcitationRatio
protected

Definition at line 114 of file G4S1Light.hh.

G4bool G4S1Light::fAlpha
protected

Definition at line 111 of file G4S1Light.hh.

G4bool G4S1Light::fExcitedNucleus
protected

Definition at line 111 of file G4S1Light.hh.

G4double G4S1Light::fKr83m
protected

Definition at line 112 of file G4S1Light.hh.

G4bool G4S1Light::fMultipleScattering
protected

Definition at line 111 of file G4S1Light.hh.

G4bool G4S1Light::fTrackSecondariesFirst
protected

Definition at line 109 of file G4S1Light.hh.

G4bool G4S1Light::fVeryHighEnergy
protected

Definition at line 111 of file G4S1Light.hh.

G4double G4S1Light::YieldFactor
protected

Definition at line 113 of file G4S1Light.hh.

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