Macros | Functions | Variables
G4S1Light.cc File Reference
#include "G4ParticleTypes.hh"
#include "G4EmProcessSubType.hh"
#include "G4Version.hh"
#include <G4SystemOfUnits.hh>
#include <G4PhysicalConstants.hh>
#include "G4S1Light.hh"

Go to the source code of this file.

Macros

#define MIN_ENE   -1*CLHEP::eV
 
#define MAX_ENE   1.*CLHEP::TeV
 
#define HIENLIM   5*CLHEP::MeV
 
#define R_TOL   0.2*CLHEP::mm
 
#define R_MAX   1*CLHEP::km
 
#define Density_LXe   2.9
 
#define Density_LAr   1.393
 
#define Density_LNe   1.207
 
#define Density_LKr   2.413
 
#define DETSIM
 

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)
 
G4double GetLiquidElectronDriftSpeed (G4double tempinput, G4double efieldinput, G4bool Miller, G4int Z)
 

Variables

G4bool diffusion = true
 
G4bool SinglePhase =false
 
G4bool ThomasImelTail =true
 
G4bool OutElectrons =true
 
G4double biExc = 0.77
 

Macro Definition Documentation

#define Density_LAr   1.393

Definition at line 78 of file G4S1Light.cc.

#define Density_LKr   2.413

Definition at line 80 of file G4S1Light.cc.

#define Density_LNe   1.207

Definition at line 79 of file G4S1Light.cc.

#define Density_LXe   2.9

Definition at line 77 of file G4S1Light.cc.

#define DETSIM
#define HIENLIM   5*CLHEP::MeV

Definition at line 59 of file G4S1Light.cc.

#define MAX_ENE   1.*CLHEP::TeV

Definition at line 58 of file G4S1Light.cc.

#define MIN_ENE   -1*CLHEP::eV

Definition at line 57 of file G4S1Light.cc.

#define R_MAX   1*CLHEP::km

Definition at line 62 of file G4S1Light.cc.

#define R_TOL   0.2*CLHEP::mm

Definition at line 61 of file G4S1Light.cc.

Function Documentation

G4int BinomFluct ( G4int  N0,
G4double  prob 
)

Definition at line 1245 of file G4S1Light.cc.

1245  {
1246  G4double mean = N0*prob;
1247  G4double sigma = sqrt(N0*prob*(1-prob));
1248  G4int N1 = 0;
1249  if ( prob == 0.00 ) return N1;
1250  if ( prob == 1.00 ) return N0;
1251 
1252  if ( N0 < 10 ) {
1253  for(G4int i = 0; i < N0; i++) {
1254  if(G4UniformRand() < prob) N1++;
1255  }
1256  }
1257  else {
1258  N1 = G4int(floor(G4RandGauss::shoot(mean,sigma)+0.5));
1259  }
1260  if ( N1 > N0 ) N1 = N0;
1261  if ( N1 < 0 ) N1 = 0;
1262  return N1;
1263 }
N1
Definition: group.cpp:53
cet::sqlite::mean
double mean(sqlite3 *db, std::string const &table_name, std::string const &column_name)
Definition: statistics.cc:16
G4double CalculateElectronLET ( G4double  E,
G4int  Z 
)

Definition at line 1218 of file G4S1Light.cc.

1218  {
1219  G4double LET;
1220  switch ( Z ) {
1221  case 54:
1222  //use a spline fit to online ESTAR data
1223  if ( E >= 1 ) LET = 58.482-61.183*log10(E)+19.749*pow(log10(E),2)+
1224  2.3101*pow(log10(E),3)-3.3469*pow(log10(E),4)+
1225  0.96788*pow(log10(E),5)-0.12619*pow(log10(E),6)+0.0065108*pow(log10(E),7);
1226  //at energies <1 keV, use a different spline, determined manually by
1227  //generating sub-keV electrons in Geant4 and looking at their ranges, since
1228  //ESTAR does not go this low
1229  else if ( E>0 && E<1 ) LET = 6.9463+815.98*E-4828*pow(E,2)+17079*pow(E,3)-
1230  36394*pow(E,4)+44553*pow(E,5)-28659*pow(E,6)+7483.8*pow(E,7);
1231  else
1232  LET = 0;
1233  break;
1234  case 18: default:
1235  if ( E >= 1 ) LET = 116.70-162.97*log10(E)+99.361*pow(log10(E),2)-
1236  33.405*pow(log10(E),3)+6.5069*pow(log10(E),4)-
1237  0.69334*pow(log10(E),5)+.031563*pow(log10(E),6);
1238  else if ( E>0 && E<1 ) LET = 100;
1239  else
1240  LET = 0;
1241  }
1242  return LET;
1243 }
cet::pow
constexpr T pow(T x)
Definition: pow.h:72
G4double GetGasElectronDriftSpeed ( G4double  efieldinput,
G4double  density 
)

Definition at line 315 of file G4S2Light.cc.

316 {
317  //Gas equation one coefficients(E/N of 1.2E-19 to 3.5E-19)
318  double gas1a=395.50266631436,gas1b=-357384143.004642,gas1c=0.518110447340587;
319  //Gas equation two coefficients(E/N of 3.5E-19 to 3.8E-17)
320  double gas2a=-592981.611357632,gas2b=-90261.9643716643,
321  gas2c=-4911.83213989609,gas2d=-115.157545835228, gas2f=-0.990440443390298,
322  gas2g=1008.30998933704,gas2h=223.711221224885;
323 
324  G4double edrift=0, gasdep=efieldinput/density, gas1fix=0, gas2fix=0;
325 
326  if ( gasdep < 1.2e-19 && gasdep >= 0 ) edrift = 4e22*gasdep;
327  if ( gasdep < 3.5e-19 && gasdep >= 1.2e-19 ) {
328  gas1fix = gas1b*pow(gasdep,gas1c); edrift = gas1a*pow(gasdep,gas1fix);
329  }
330  if ( gasdep < 3.8e-17 && gasdep >= 3.5e-19 ) {
331  gas2fix = log(gas2g*gasdep);
332  edrift = (gas2a+gas2b*gas2fix+gas2c*pow(gas2fix,2)+gas2d*pow(gas2fix,3)+
333  gas2f*pow(gas2fix,4))*(gas2h*exp(gasdep));
334  }
335  if ( gasdep >= 3.8e-17 ) edrift = 6e21*gasdep-32279;
336 
337  return 0.5*EMASS*pow(edrift*(CLHEP::cm/CLHEP::s),2.);
338 }
genie::units::cm
static constexpr double cm
Definition: Units.h:68
EMASS
#define EMASS
Definition: G4S1Light.hh:25
cet::pow
constexpr T pow(T x)
Definition: pow.h:72
e
const double e
Definition: gUpMuFluxGen.cxx:165
s
static QCString * s
Definition: config.cpp:1042
G4double GetLiquidElectronDriftSpeed ( double  T,
double  F,
G4bool  M,
G4int  Z 
)
G4double GetLiquidElectronDriftSpeed ( G4double  tempinput,
G4double  efieldinput,
G4bool  Miller,
G4int  Z 
)

Definition at line 1118 of file G4S1Light.cc.

1119  {
1120  if(efieldinput<0) efieldinput *= (-1);
1121  //Liquid equation one (165K) coefficients
1122  G4double onea=144623.235704015,
1123  oneb=850.812714257629,
1124  onec=1192.87056676815,
1125  oned=-395969.575204061,
1126  onef=-355.484170008875,
1127  oneg=-227.266219627672,
1128  oneh=223831.601257495,
1129  onei=6.1778950907965,
1130  onej=18.7831533426398,
1131  onek=-76132.6018884368;
1132  //Liquid equation two (200K) coefficients
1133  G4double twoa=17486639.7118995,
1134  twob=-113.174284723134,
1135  twoc=28.005913193763,
1136  twod=167994210.094027,
1137  twof=-6766.42962575088,
1138  twog=901.474643115395,
1139  twoh=-185240292.471665,
1140  twoi=-633.297790813084,
1141  twoj=87.1756135457949;
1142  //Liquid equation three (230K) coefficients
1143  G4double thra=10626463726.9833,
1144  thrb=224025158.134792,
1145  thrc=123254826.300172,
1146  thrd=-4563.5678061122,
1147  thrf=-1715.269592063,
1148  thrg=-694181.921834368,
1149  thrh=-50.9753281079838,
1150  thri=58.3785811395493,
1151  thrj=201512.080026704;
1152  G4double y1=0,y2=0,f1=0,f2=0,f3=0,edrift=0,
1153  t1=0,t2=0,frac=0,slope=0,intercept=0;
1154 
1155  //Equations defined
1156  f1=onea/(1+exp(-(efieldinput-oneb)/onec))+oned/
1157  (1+exp(-(efieldinput-onef)/oneg))+
1158  oneh/(1+exp(-(efieldinput-onei)/onej))+onek;
1159  f2=twoa/(1+exp(-(efieldinput-twob)/twoc))+twod/
1160  (1+exp(-(efieldinput-twof)/twog))+
1161  twoh/(1+exp(-(efieldinput-twoi)/twoj));
1162  f3=thra*exp(-thrb*efieldinput)+thrc*exp(-(pow(efieldinput-thrd,2))/
1163  (thrf*thrf))+
1164  thrg*exp(-(pow(efieldinput-thrh,2)/(thri*thri)))+thrj;
1165 
1166  if(efieldinput<20 && efieldinput>=0) {
1167  f1=2951*efieldinput;
1168  f2=5312*efieldinput;
1169  f3=7101*efieldinput;
1170  }
1171  //Cases for tempinput decides which 2 equations to use lin. interpolation
1172  if(tempinput<200.0 && tempinput>165.0) {
1173  y1=f1;
1174  y2=f2;
1175  t1=165.0;
1176  t2=200.0;
1177  }
1178  if(tempinput<230.0 && tempinput>200.0) {
1179  y1=f2;
1180  y2=f3;
1181  t1=200.0;
1182  t2=230.0;
1183  }
1184  if((tempinput>230.0 || tempinput<165.0) && !Miller) {
1185  G4cout << "\nWARNING: TEMPERATURE OUT OF RANGE (165-230 K)\n";
1186  return 0;
1187  }
1188  if (tempinput == 165.0) edrift = f1;
1189  else if (tempinput == 200.0) edrift = f2;
1190  else if (tempinput == 230.0) edrift = f3;
1191  else { //Linear interpolation
1192  frac=(tempinput-t1)/(t2-t1);
1193  slope = (y1-y2)/(t1-t2);
1194  intercept=y1-slope*t1;
1195  edrift=slope*tempinput+intercept;
1196  }
1197 
1198  if ( Miller ) {
1199  if ( efieldinput <= 40. )
1200  edrift = -0.13274+0.041082*efieldinput-0.0006886*pow(efieldinput,2.)+
1201  5.5503e-6*pow(efieldinput,3.);
1202  else
1203  edrift = 0.060774*efieldinput/pow(1+0.11336*pow(efieldinput,0.5218),2.);
1204  if ( efieldinput >= 1e5 ) edrift = 2.7;
1205  if ( efieldinput >= 100 )
1206  edrift -= 0.017 * ( tempinput - 163 );
1207  else
1208  edrift += 0.017 * ( tempinput - 163 );
1209  edrift *= 1e5; //put into units of cm/sec. from mm/usec.
1210  }
1211  if ( Z == 18 ) edrift = 1e5 * (
1212  .097384*pow(log10(efieldinput),3.0622)-.018614*sqrt(efieldinput) );
1213  if ( edrift < 0 ) edrift = 0.;
1214  edrift = 0.5*EMASS*pow(edrift*CLHEP::cm/CLHEP::s,2.);
1215  return edrift;
1216 }
genie::units::cm
static constexpr double cm
Definition: Units.h:68
EMASS
#define EMASS
Definition: G4S1Light.hh:25
reco_momentum_tuples.f2
f2
Definition: reco_momentum_tuples.py:79
cet::pow
constexpr T pow(T x)
Definition: pow.h:72
test_gpu_visible.t1
t1
Definition: test_gpu_visible.py:16
reco_momentum_tuples.f3
f3
Definition: reco_momentum_tuples.py:80
s
static QCString * s
Definition: config.cpp:1042
reco_momentum_tuples.f1
f1
Definition: reco_momentum_tuples.py:78
void InitMatPropValues ( G4MaterialPropertiesTable *  nobleElementMat)

Definition at line 1265 of file G4S1Light.cc.

1265  {
1266  char xCoord[80]; char yCoord[80]; char zCoord[80];
1267  char numExc[80]; char numIon[80]; char numPho[80]; char numEle[80];
1268  char trackL[80]; char time00[80]; char time01[80]; char energy[80];
1269 
1270  // for loop to initialize the interaction site mat'l properties
1271  for( G4int i=0; i<10000; i++ ) {
1272  sprintf(xCoord,"POS_X_%d",i); sprintf(yCoord,"POS_Y_%d",i);
1273  sprintf(zCoord,"POS_Z_%d",i);
1274  nobleElementMat->AddConstProperty( xCoord, 999*CLHEP::km );
1275  nobleElementMat->AddConstProperty( yCoord, 999*CLHEP::km );
1276  nobleElementMat->AddConstProperty( zCoord, 999*CLHEP::km );
1277  sprintf(numExc,"N_EXC_%d",i); sprintf(numIon,"N_ION_%d",i);
1278  sprintf(numPho,"N_PHO_%d",i); sprintf(numEle,"N_ELE_%d",i);
1279  nobleElementMat->AddConstProperty( numExc, 0 );
1280  nobleElementMat->AddConstProperty( numIon, 0 );
1281  nobleElementMat->AddConstProperty( numPho, 0 );
1282  nobleElementMat->AddConstProperty( numEle, 0 );
1283  sprintf(trackL,"TRACK_%d",i); sprintf(energy,"ENRGY_%d",i);
1284  sprintf(time00,"TIME0_%d",i); sprintf(time01,"TIME1_%d",i);
1285  nobleElementMat->AddConstProperty( trackL, 0*CLHEP::um );
1286  nobleElementMat->AddConstProperty( energy, 0*CLHEP::eV );
1287  nobleElementMat->AddConstProperty( time00, DBL_MAX );
1288  nobleElementMat->AddConstProperty( time01,-1*CLHEP::ns );
1289  }
1290 
1291  // we initialize the total number of interaction sites, a variable for
1292  // updating the amount of energy deposited thus far in the medium, and a
1293  // variable for storing the amount of energy expected to be deposited
1294  nobleElementMat->AddConstProperty( "TOTALNUM_INT_SITES", 0 );
1295  nobleElementMat->AddConstProperty( "ENERGY_DEPOSIT_TOT", 0*CLHEP::keV );
1296  nobleElementMat->AddConstProperty( "ENERGY_DEPOSIT_GOL", 0*CLHEP::MeV );
1297  return;
1298 }
genie::units::keV
static constexpr double keV
Definition: Units.h:128
generate_CCQE_events.energy
string energy
Definition: generate_CCQE_events.py:15
genie::units::MeV
static constexpr double MeV
Definition: Units.h:129
genie::units::km
static constexpr double km
Definition: Units.h:64
genie::units::eV
static constexpr double eV
Definition: Units.h:127
ns
QAsciiDict< Entry > ns
Definition: tclscanner.cpp:1147
G4double UnivScreenFunc ( G4double  E,
G4double  Z,
G4double  A 
)

Definition at line 1300 of file G4S1Light.cc.

1300  {
1301  G4double a_0 = 5.29e-11*CLHEP::m; G4double a = 0.626*a_0*pow(Z,(-1./3.));
1302  G4double epsilon_0 = 8.854e-12*(CLHEP::farad/CLHEP::m);
1303  G4double epsilon = (a*E*2.*CLHEP::twopi*epsilon_0)/(2*pow(CLHEP::eplus,2.)*pow(Z,2.));
1304  G4double zeta_0 = pow(Z,(1./6.)); G4double m_N = A*1.66e-27*CLHEP::kg;
1305  G4double hbar = 6.582e-16*CLHEP::eV*CLHEP::s;
1306  if ( Z == 54 ) {
1307  epsilon *= 1.068; //zeta_0 = 1.63;
1308  } //special case for LXe from Bezrukov et al. 2011
1309  G4double s_n = log(1+1.1383*epsilon)/(2.*(epsilon +
1310  0.01321*pow(epsilon,0.21226) +
1311  0.19593*sqrt(epsilon)));
1312  G4double s_e = (a_0*zeta_0/a)*hbar*sqrt(8*epsilon*2.*CLHEP::twopi*epsilon_0/
1313  (a*m_N*pow(CLHEP::eplus,2.)));
1314  return 1.38e5*0.5*(1+tanh(50*epsilon-0.25))*epsilon*(s_e/s_n);
1315 }
hbar
const double hbar
Definition: CommonIMapIncludes.h:25
cet::pow
constexpr T pow(T x)
Definition: pow.h:72
genie::units::kg
static constexpr double kg
Definition: Units.h:143
genie::units::eV
static constexpr double eV
Definition: Units.h:127
a
const double a
Definition: gUpMuFluxGen.cxx:164
wirecell.units.twopi
int twopi
Definition: units.py:12
cache_state.m
m
Definition: cache_state.py:415
s
static QCString * s
Definition: config.cpp:1042

Variable Documentation

G4double biExc = 0.77

Definition at line 70 of file G4S1Light.cc.

G4bool diffusion = true

Definition at line 63 of file G4S1Light.cc.

G4bool OutElectrons =true

Definition at line 65 of file G4S1Light.cc.

G4bool SinglePhase =false

Definition at line 65 of file G4S1Light.cc.

G4bool ThomasImelTail =true

Definition at line 65 of file G4S1Light.cc.

Generated by   doxygen 1.8.11