All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | Private Member Functions | Private Attributes | List of all members
NestAlg Class Reference

#include <NestAlg.h>

Public Member Functions

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

Private Member Functions

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

Private Attributes

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

Detailed Description

Definition at line 14 of file NestAlg.h.

Constructor & Destructor Documentation

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

Definition at line 49 of file NestAlg.cxx.

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

Definition at line 65 of file NestAlg.cxx.

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

Member Function Documentation

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

Definition at line 1189 of file NestAlg.cxx.

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

Definition at line 1161 of file NestAlg.cxx.

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

Definition at line 81 of file NestAlg.cxx.

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

Definition at line 28 of file NestAlg.h.

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

Definition at line 1048 of file NestAlg.cxx.

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

Definition at line 1059 of file NestAlg.cxx.

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

Definition at line 1213 of file NestAlg.cxx.

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

Definition at line 27 of file NestAlg.h.

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

Definition at line 26 of file NestAlg.h.

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

Definition at line 25 of file NestAlg.h.

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

Definition at line 24 of file NestAlg.h.

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

Member Data Documentation

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

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

Definition at line 46 of file NestAlg.h.

double NestAlg::fEnergyDep
private

energy deposited by the step

Definition at line 44 of file NestAlg.h.

CLHEP::HepRandomEngine& NestAlg::fEngine
private

random engine

Definition at line 49 of file NestAlg.h.

double NestAlg::fExcitationRatio
private

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

Definition at line 40 of file NestAlg.h.

int NestAlg::fNumIonElectrons
private

number of ionization electrons produced by step

Definition at line 43 of file NestAlg.h.

int NestAlg::fNumScintPhotons
private

number of photons produced by the step

Definition at line 42 of file NestAlg.h.

G4VParticleChange NestAlg::fParticleChange
private

pointer to G4VParticleChange

Definition at line 45 of file NestAlg.h.

double 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.5