1 #define AVO 6.022e23 //Avogadro's number (#/mol) 2 #define EMASS 9.109e-31 * CLHEP::kg 3 #define MillerDriftSpeed true 5 #define GASGAP 0.25 * CLHEP::cm //S2 generation region 6 #define BORDER 0 * CLHEP::cm //liquid-gas border z-coordinate 8 #define QE_EFF 1 //a base or maximum quantum efficiency 9 #define phe_per_e 1 //S2 gain for quick studies 12 #define WIN 0 * CLHEP::mm //top Cu block (also, quartz window) 13 #define TOP 0 //top grid wires 14 #define ANE 0 //anode mesh 15 #define SRF 0 //liquid-gas interface 16 #define GAT 0 //gate grid 17 #define CTH 0 //cathode grid 18 #define BOT 0 //bottom PMT grid 19 #define PMT 0 //bottom Cu block and PMTs 20 #define MIN_ENE -1 * CLHEP::eV //lets you turn NEST off BELOW a certain energy 21 #define MAX_ENE 1. * CLHEP::TeV //lets you turn NEST off ABOVE a certain energy 22 #define HIENLIM 5 * CLHEP::MeV //energy at which Doke model used exclusively 23 #define R_TOL 0.2 * CLHEP::mm //tolerance (for edge events) 24 #define R_MAX 1 * CLHEP::km //for corraling diffusing electrons 25 #define Density_LXe 2.9 //reference density for density-dep. effects 26 #define Density_LAr 1.393 27 #define Density_LNe 1.207 28 #define Density_LKr 2.413 30 #include "Geant4/G4Ions.hh" 31 #include "Geant4/G4OpticalPhoton.hh" 32 #include "Geant4/G4VProcess.hh" 37 #include "CLHEP/Random/RandFlat.h" 38 #include "CLHEP/Random/RandGauss.h" 84 CLHEP::RandGauss GaussGen(
fEngine);
85 CLHEP::RandFlat UniformGen(
fEngine);
97 bool fTrackSecondariesFirst =
false;
98 bool fExcitedNucleus =
false;
99 bool fVeryHighEnergy =
false;
101 bool fMultipleScattering =
false;
104 if (aTrack.GetParentID() == 0 && aTrack.GetCurrentStepNumber() == 1) {
105 fExcitedNucleus =
false;
106 fVeryHighEnergy =
false;
108 fMultipleScattering =
false;
111 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
112 G4ParticleDefinition* pDef = aParticle->GetDefinition();
113 G4String particleName = pDef->GetParticleName();
114 const G4Material* aMaterial = aStep.GetPreStepPoint()->GetMaterial();
115 const G4Material* bMaterial = aStep.GetPostStepPoint()->GetMaterial();
117 if ((particleName ==
"neutron" || particleName ==
"antineutron") &&
118 aStep.GetTotalEnergyDeposit() <= 0)
126 G4Element *ElementA = NULL, *ElementB = NULL;
128 const G4ElementVector* theElementVector1 = aMaterial->GetElementVector();
129 ElementA = (*theElementVector1)[0];
132 const G4ElementVector* theElementVector2 = bMaterial->GetElementVector();
133 ElementB = (*theElementVector2)[0];
136 G4bool NobleNow =
false, NobleLater =
false;
138 z1 = (G4int)(ElementA->GetZ());
142 z2 = (G4int)(ElementB->GetZ());
145 if (z1 == 2 || z1 == 10 || z1 == 18 || z1 == 36 || z1 == 54) {
155 if (z2 == 2 || z2 == 10 || z2 == 18 || z2 == 36 || z2 == 54) {
170 G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
171 G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
172 G4ThreeVector
x1 = pPostStepPoint->GetPosition();
173 G4ThreeVector x0 = pPreStepPoint->GetPosition();
174 G4double evtStrt = pPreStepPoint->GetGlobalTime();
175 G4double
t0 = pPreStepPoint->GetLocalTime();
176 G4double
t1 = pPostStepPoint->GetLocalTime();
182 G4bool outside =
false, inside =
false, InsAndOuts =
false;
183 G4MaterialPropertiesTable* aMaterialPropertiesTable = aMaterial->GetMaterialPropertiesTable();
184 if (NobleNow && !NobleLater) outside =
true;
185 if (!NobleNow && NobleLater) {
186 aMaterial = bMaterial;
189 aMaterialPropertiesTable = bMaterial->GetMaterialPropertiesTable();
191 if (NobleNow && NobleLater && aMaterial->GetDensity() != bMaterial->GetDensity())
195 G4double Density = aMaterial->GetDensity() / (CLHEP::g / CLHEP::cm3);
196 G4double nDensity = Density *
AVO;
197 G4int Phase = aMaterial->GetState();
198 G4double ElectricField(0.), FieldSign(0.);
199 G4bool GlobalFields =
false;
202 ElectricField = aMaterialPropertiesTable->GetConstProperty(
"ELECTRICFIELD");
206 if (x1[2] <
WIN && x1[2] >
TOP && Phase == kStateGas)
207 ElectricField = aMaterialPropertiesTable->GetConstProperty(
"ELECTRICFIELDWINDOW");
208 else if (x1[2] <
TOP && x1[2] >
ANE && Phase == kStateGas)
209 ElectricField = aMaterialPropertiesTable->GetConstProperty(
"ELECTRICFIELDTOP");
210 else if (x1[2] <
ANE && x1[2] >
SRF && Phase == kStateGas)
211 ElectricField = aMaterialPropertiesTable->GetConstProperty(
"ELECTRICFIELDANODE");
212 else if (x1[2] <
SRF && x1[2] >
GAT && Phase == kStateLiquid)
213 ElectricField = aMaterialPropertiesTable->GetConstProperty(
"ELECTRICFIELDSURFACE");
214 else if (x1[2] <
GAT && x1[2] >
CTH && Phase == kStateLiquid)
215 ElectricField = aMaterialPropertiesTable->GetConstProperty(
"ELECTRICFIELDGATE");
216 else if (x1[2] <
CTH && x1[2] >
BOT && Phase == kStateLiquid)
217 ElectricField = aMaterialPropertiesTable->GetConstProperty(
"ELECTRICFIELDCATHODE");
218 else if (x1[2] <
BOT && x1[2] >
PMT && Phase == kStateLiquid)
219 ElectricField = aMaterialPropertiesTable->GetConstProperty(
"ELECTRICFIELDBOTTOM");
221 ElectricField = aMaterialPropertiesTable->GetConstProperty(
"ELECTRICFIELD");
223 if (ElectricField >= 0)
228 G4double Temperature = aMaterial->GetTemperature();
229 G4double ScintillationYield, ResolutionScale, R0 = 1.0 * CLHEP::um, DokeBirks[3],
230 ThomasImel = 0.00, delta = 1 * CLHEP::mm;
233 G4double PhotMean = 7 * CLHEP::eV, PhotWidth = 1.0 * CLHEP::eV;
234 G4double SingTripRatioR, SingTripRatioX, tau1, tau3, tauR = 0 * CLHEP::ns;
237 ScintillationYield = 1 / (41.3 * CLHEP::eV);
239 ResolutionScale = 0.2;
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);
246 ScintillationYield = 1 / (29.2 * CLHEP::eV);
248 ResolutionScale = 0.13;
249 PhotMean = 15.5 * CLHEP::eV;
250 PhotWidth = 0.26 * CLHEP::eV;
251 tau1 = GaussGen.fire(10.0 * CLHEP::ns, 10. * CLHEP::ns);
252 tau3 = GaussGen.fire(15.4e3 * CLHEP::ns, 200 * CLHEP::ns);
255 ScintillationYield = 1 / (19.5 * CLHEP::eV);
257 ResolutionScale = 0.107;
258 R0 = 1.568 * CLHEP::um;
260 ThomasImel = 0.156977 * pow(ElectricField, -0.1);
261 DokeBirks[0] = 0.07 * pow((ElectricField / 1.0e3), -0.85);
266 DokeBirks[0] = 0.0003;
270 PhotMean = 9.69 * CLHEP::eV;
271 PhotWidth = 0.22 * CLHEP::eV;
272 tau1 = GaussGen.fire(6.5 * CLHEP::ns, 0.8 * CLHEP::ns);
273 tau3 = GaussGen.fire(1300 * CLHEP::ns, 50 * CLHEP::ns);
274 tauR = GaussGen.fire(0.8 * CLHEP::ns, 0.2 * CLHEP::ns);
278 if (Phase == kStateGas)
279 ScintillationYield = 1 / (30.0 * CLHEP::eV);
281 ScintillationYield = 1 / (15.0 * CLHEP::eV);
283 ResolutionScale = 0.05;
284 PhotMean = 8.43 * CLHEP::eV;
285 tau1 = GaussGen.fire(2.8 * CLHEP::ns, .04 * CLHEP::ns);
286 tau3 = GaussGen.fire(93. * CLHEP::ns, 1.1 * CLHEP::ns);
287 tauR = GaussGen.fire(12. * CLHEP::ns, .76 * CLHEP::ns);
292 ScintillationYield = 48.814 + 0.80877 * Density + 2.6721 * pow(Density, 2.);
293 ScintillationYield /= CLHEP::keV;
297 ResolutionScale = 1.00 *
298 (0.12724 - 0.032152 * Density - 0.0013492 * pow(Density, 2.));
300 if (Phase == kStateLiquid) {
301 ResolutionScale *= 1.5;
302 R0 = 16.6 * CLHEP::um;
305 R0 = 69.492 * pow(ElectricField, -0.50422) * CLHEP::um;
307 DokeBirks[0] = 19.171 * pow(ElectricField + 25.552, -0.83057) + 0.026772;
317 delta = 0.4 * CLHEP::mm;
318 PhotMean = 6.97 * CLHEP::eV;
319 PhotWidth = 0.23 * CLHEP::eV;
324 tau1 = GaussGen.fire(3.1 * CLHEP::ns, .7 * CLHEP::ns);
325 tau3 = GaussGen.fire(24. * CLHEP::ns, 1. * CLHEP::ns);
327 else if (Phase == kStateGas) {
332 ScintillationYield = 1 / (12.98 * CLHEP::eV);
334 R0 = 0.0 * CLHEP::um;
335 G4double Townsend = (ElectricField / nDensity) * 1e17;
336 DokeBirks[0] = 0.0000;
337 DokeBirks[2] = 0.1933 * pow(Density, 2.6199) + 0.29754 -
338 (0.045439 * pow(Density, 2.4689) + 0.066034) * log10(ElectricField);
339 if (ElectricField > 6990) DokeBirks[2] = 0.0;
340 if (ElectricField < 1000) DokeBirks[2] = 0.2;
341 if (ElectricField < 100.) {
346 ThomasImel = 0.041973 * pow(Density, 1.8105);
347 else if (Density >= 0.061 && Density <= 0.167)
348 ThomasImel = 5.9583e-5 + 0.0048523 * Density - 0.023109 * pow(Density, 2.);
350 ThomasImel = 6.2552e-6 * pow(Density, -1.9963);
351 if (ElectricField) ThomasImel = 1.2733e-5 * pow(Townsend / Density, -0.68426);
353 PhotMean = 7.1 * CLHEP::eV;
354 PhotWidth = 0.2 * CLHEP::eV;
355 tau1 = GaussGen.fire(5.18 * CLHEP::ns, 1.55 * CLHEP::ns);
356 tau3 = GaussGen.fire(100.1 * CLHEP::ns, 7.9 * CLHEP::ns);
359 tau1 = 3.5 * CLHEP::ns;
360 tau3 = 20. * CLHEP::ns;
361 tauR = 40. * CLHEP::ns;
366 G4double anExcitationEnergy =
367 ((
const G4Ions*)(pDef))->GetExcitationEnergy();
368 G4double TotalEnergyDeposit =
369 aMaterialPropertiesTable->GetConstProperty(
"ENERGY_DEPOSIT_TOT");
370 G4bool convert =
false, annihil =
false;
372 if (pPreStepPoint->GetKineticEnergy() >= (2 * CLHEP::electron_mass_c2) &&
373 !pPostStepPoint->GetKineticEnergy() && !aStep.GetTotalEnergyDeposit() &&
374 aParticle->GetPDGcode() == 22) {
376 TotalEnergyDeposit = CLHEP::electron_mass_c2;
378 if (pPreStepPoint->GetKineticEnergy() && !pPostStepPoint->GetKineticEnergy() &&
379 aParticle->GetPDGcode() == -11) {
381 TotalEnergyDeposit += aStep.GetTotalEnergyDeposit();
383 G4bool either =
false;
384 if (inside || outside || convert || annihil || InsAndOuts) either =
true;
386 if (anExcitationEnergy < 100 * CLHEP::eV && aStep.GetTotalEnergyDeposit() < 1 * CLHEP::eV &&
387 !either && !fExcitedNucleus)
390 if (!annihil) TotalEnergyDeposit += aStep.GetTotalEnergyDeposit();
392 aMaterialPropertiesTable->AddConstProperty(
"ENERGY_DEPOSIT_TOT", TotalEnergyDeposit);
394 TotalEnergyDeposit = aStep.GetTotalEnergyDeposit();
399 G4double InitialKinetEnergy = aMaterialPropertiesTable->GetConstProperty(
"ENERGY_DEPOSIT_GOL");
401 if (InitialKinetEnergy == 0) {
402 G4double tE = pPreStepPoint->GetKineticEnergy() + anExcitationEnergy;
403 if ((fabs(tE - 1.8 * CLHEP::keV) < CLHEP::eV || fabs(tE - 9.4 * CLHEP::keV) < CLHEP::eV) &&
404 Phase == kStateLiquid && z1 == 54)
405 tE = 9.4 * CLHEP::keV;
406 if (fKr83m && ElectricField != 0) DokeBirks[2] = 0.10;
407 aMaterialPropertiesTable->AddConstProperty(
"ENERGY_DEPOSIT_GOL", tE);
411 if (anExcitationEnergy) fExcitedNucleus =
true;
416 aMaterialPropertiesTable->AddConstProperty(
417 "ENERGY_DEPOSIT_GOL", InitialKinetEnergy - pPostStepPoint->GetKineticEnergy());
418 if (aMaterialPropertiesTable->GetConstProperty(
"ENERGY_DEPOSIT_GOL") < 0)
419 aMaterialPropertiesTable->AddConstProperty(
"ENERGY_DEPOSIT_GOL", 0);
424 aMaterialPropertiesTable->AddConstProperty(
425 "ENERGY_DEPOSIT_GOL", InitialKinetEnergy + pPreStepPoint->GetKineticEnergy());
426 if (TotalEnergyDeposit > 0 && InitialKinetEnergy == 0) {
427 aMaterialPropertiesTable->AddConstProperty(
"ENERGY_DEPOSIT_GOL", 0);
428 TotalEnergyDeposit = .000000;
434 aMaterialPropertiesTable->AddConstProperty(
"ENERGY_DEPOSIT_GOL",
435 (-0.1 * CLHEP::keV) + InitialKinetEnergy -
436 pPostStepPoint->GetKineticEnergy());
438 bMaterial->GetMaterialPropertiesTable()->GetConstProperty(
"ENERGY_DEPOSIT_GOL");
439 bMaterial->GetMaterialPropertiesTable()->AddConstProperty(
440 "ENERGY_DEPOSIT_GOL",
441 (-0.1 * CLHEP::keV) + InitialKinetEnergy + pPreStepPoint->GetKineticEnergy());
442 if (aMaterialPropertiesTable->GetConstProperty(
"ENERGY_DEPOSIT_GOL") < 0)
443 aMaterialPropertiesTable->AddConstProperty(
"ENERGY_DEPOSIT_GOL", 0);
444 if (bMaterial->GetMaterialPropertiesTable()->GetConstProperty(
"ENERGY_DEPOSIT_GOL") < 0)
445 bMaterial->GetMaterialPropertiesTable()->AddConstProperty(
"ENERGY_DEPOSIT_GOL", 0);
448 aMaterialPropertiesTable->GetConstProperty(
"ENERGY_DEPOSIT_GOL");
450 InitialKinetEnergy += 2 * CLHEP::electron_mass_c2;
453 if (convert) InitialKinetEnergy -= 2 * CLHEP::electron_mass_c2;
455 aMaterialPropertiesTable->AddConstProperty(
"ENERGY_DEPOSIT_GOL", InitialKinetEnergy);
456 if (anExcitationEnergy < 1
e-100 && aStep.GetTotalEnergyDeposit() == 0 &&
457 aMaterialPropertiesTable->GetConstProperty(
"ENERGY_DEPOSIT_GOL") == 0 &&
458 aMaterialPropertiesTable->GetConstProperty(
"ENERGY_DEPOSIT_TOT") == 0)
462 if (aTrack.GetCreatorProcess())
463 procName = aTrack.GetCreatorProcess()->GetProcessName();
466 if (procName ==
"eBrem" && outside && !
OutElectrons) fMultipleScattering =
true;
469 if (fAlpha) delta = 1000. * CLHEP::km;
473 if (aParticle->GetPDGcode() == 11 && !
OutElectrons) fMultipleScattering =
true;
480 G4bool exists =
false;
481 for (i = 0; i < j; i++) {
483 sprintf(xCoord,
"POS_X_%d", i);
484 sprintf(yCoord,
"POS_Y_%d", i);
485 sprintf(zCoord,
"POS_Z_%d", i);
486 pos[0] = aMaterialPropertiesTable->GetConstProperty(xCoord);
487 pos[1] = aMaterialPropertiesTable->GetConstProperty(yCoord);
488 pos[2] = aMaterialPropertiesTable->GetConstProperty(zCoord);
489 if (sqrt(pow(x1[0] - pos[0], 2.) + pow(x1[1] - pos[1], 2.) + pow(x1[2] - pos[2], 2.)) < delta) {
494 if (!exists && TotalEnergyDeposit) {
496 sprintf(xCoord,
"POS_X_%i", j);
497 sprintf(yCoord,
"POS_Y_%i", j);
498 sprintf(zCoord,
"POS_Z_%i", j);
500 aMaterialPropertiesTable->AddConstProperty(xCoord, x1[0]);
501 aMaterialPropertiesTable->AddConstProperty(yCoord, x1[1]);
502 aMaterialPropertiesTable->AddConstProperty(zCoord, x1[2]);
504 aMaterialPropertiesTable->
505 AddConstProperty(
"TOTALNUM_INT_SITES", j);
513 G4double
a1 = ElementA->GetA();
514 z2 = pDef->GetAtomicNumber();
515 G4double
a2 = (G4double)(pDef->GetAtomicMass());
516 if (particleName ==
"alpha" || (z2 == 2 && a2 == 4))
518 if (fAlpha ||
abs(aParticle->GetPDGcode()) == 2112) a2 =
a1;
519 G4double epsilon = 11.5 * (TotalEnergyDeposit / CLHEP::keV) * pow(z1, (-7. / 3.));
520 G4double gamma = 3. * pow(epsilon, 0.15) + 0.7 * pow(epsilon, 0.6) + epsilon;
521 G4double kappa = 0.133 * pow(z1, (2. / 3.)) * pow(a2, (-1. / 2.)) * (2. / 3.);
523 if ((z1 == z2 && pDef->GetParticleType() ==
"nucleus") || particleName ==
"neutron" ||
524 particleName ==
"antineutron") {
526 if (z1 == 18 && Phase == kStateLiquid)
531 if (ElectricField == 0 && Phase == kStateLiquid) {
532 if (z1 == 54) ThomasImel = 0.19;
533 if (z1 == 18) ThomasImel = 0.25;
535 fExcitationRatio = 0.69337 + 0.3065 * exp(-0.008806 * pow(ElectricField, 0.76313));
541 G4double MeanNumberOfQuanta =
542 ScintillationYield * TotalEnergyDeposit;
545 G4double sigma = sqrt(ResolutionScale * MeanNumberOfQuanta);
547 G4int(floor(GaussGen.fire(MeanNumberOfQuanta, sigma) + 0.5));
549 if (LeffVar > 1) { LeffVar = 1.00000; }
550 if (LeffVar < 0) { LeffVar = 0; }
554 if (TotalEnergyDeposit < 1 / ScintillationYield || NumQuanta < 0) NumQuanta = 0;
558 G4int NumIons = NumQuanta - NumExcitons;
564 G4double dE, dx = 0, LET = 0, recombProb;
565 dE = TotalEnergyDeposit / CLHEP::MeV;
566 if (particleName !=
"e-" && particleName !=
"e+" && z1 != z2 && particleName !=
"mu-" &&
567 particleName !=
"mu+") {
573 if (LET) dx = dE / (Density * LET);
574 if (
abs(aParticle->GetPDGcode()) == 2112) dx = 0;
577 dx = aStep.GetStepLength() / CLHEP::cm;
578 if (dx) LET = (dE / dx) * (1 / Density);
579 if (LET > 0 && dE > 0 && dx > 0) {
581 if (j == 1 && ratio < 0.7 && !
ThomasImelTail && particleName ==
"e-") {
587 DokeBirks[1] = DokeBirks[0] / (1 - DokeBirks[2]);
589 recombProb = (DokeBirks[0] * LET) / (1 + DokeBirks[1] * LET) + DokeBirks[2];
590 if (Phase == kStateLiquid) {
591 if (z1 == 54) recombProb *= (Density /
Density_LXe);
592 if (z1 == 18) recombProb *= (Density /
Density_LAr);
595 if (recombProb < 0) recombProb = 0;
596 if (recombProb > 1) recombProb = 1;
600 G4int NumPhotons = NumExcitons +
BinomFluct(NumIons, recombProb);
601 G4int NumElectrons = NumQuanta - NumPhotons;
616 sprintf(numExc,
"N_EXC_%i", counter);
617 sprintf(numIon,
"N_ION_%i", counter);
618 aMaterialPropertiesTable->AddConstProperty(numExc, NumExcitons);
619 aMaterialPropertiesTable->AddConstProperty(numIon, NumIons);
620 sprintf(numPho,
"N_PHO_%i", counter);
621 sprintf(numEle,
"N_ELE_%i", counter);
622 aMaterialPropertiesTable->AddConstProperty(numPho, NumPhotons);
623 aMaterialPropertiesTable->AddConstProperty(numEle, NumElectrons);
631 sprintf(trackL,
"TRACK_%i", counter);
632 sprintf(energy,
"ENRGY_%i", counter);
633 sprintf(time00,
"TIME0_%i", counter);
634 sprintf(time01,
"TIME1_%i", counter);
635 delta = aMaterialPropertiesTable->GetConstProperty(trackL);
636 G4double energ = aMaterialPropertiesTable->GetConstProperty(energy);
637 delta += dx * CLHEP::cm;
638 energ += dE * CLHEP::MeV;
639 aMaterialPropertiesTable->AddConstProperty(trackL, delta);
640 aMaterialPropertiesTable->AddConstProperty(energy, energ);
641 if (TotalEnergyDeposit > 0) {
642 G4double deltaTime = aMaterialPropertiesTable->GetConstProperty(time00);
645 if (aParticle->GetCharge() != 0) {
646 if (t0 < deltaTime) aMaterialPropertiesTable->AddConstProperty(time00, t0);
649 if (t1 < deltaTime) aMaterialPropertiesTable->AddConstProperty(time00, t1);
651 deltaTime = aMaterialPropertiesTable->GetConstProperty(time01);
653 if (t1 > deltaTime) aMaterialPropertiesTable->AddConstProperty(time01, t1);
658 aMaterialPropertiesTable->GetConstProperty(
"ENERGY_DEPOSIT_TOT");
660 aMaterialPropertiesTable->GetConstProperty(
"ENERGY_DEPOSIT_GOL");
661 if (InitialKinetEnergy >
HIENLIM &&
abs(aParticle->GetPDGcode()) != 2112) fVeryHighEnergy =
true;
663 if (fVeryHighEnergy && !fExcitedNucleus)
664 safety = 0.2 * CLHEP::keV;
666 safety = 2. * CLHEP::eV;
669 if (!anExcitationEnergy && pDef->GetParticleType() ==
"nucleus" &&
670 aTrack.GetTrackStatus() != fAlive && !fAlpha)
671 InitialKinetEnergy = TotalEnergyDeposit;
672 if (particleName ==
"neutron" || particleName ==
"antineutron")
673 InitialKinetEnergy = TotalEnergyDeposit;
678 if (
std::abs(TotalEnergyDeposit - InitialKinetEnergy) < safety ||
679 TotalEnergyDeposit >= InitialKinetEnergy) {
686 for (i = 0; i < j; i++) {
687 sprintf(numPho,
"N_PHO_%d", i);
688 sprintf(numEle,
"N_ELE_%d", i);
689 sprintf(trackL,
"TRACK_%d", i);
690 sprintf(energy,
"ENRGY_%d", i);
692 dx += aMaterialPropertiesTable->GetConstProperty(trackL);
693 dE += aMaterialPropertiesTable->GetConstProperty(energy);
697 if (fVeryHighEnergy) buffer = 1;
698 fParticleChange.SetNumberOfSecondaries(buffer * (NumPhotons + NumElectrons));
700 if (fTrackSecondariesFirst) {
701 if (aTrack.GetTrackStatus() == fAlive)
fParticleChange.ProposeTrackStatus(fSuspend);
705 for (i = 0; i < j; i++) {
708 sprintf(xCoord,
"POS_X_%d", i);
709 sprintf(yCoord,
"POS_Y_%d", i);
710 sprintf(zCoord,
"POS_Z_%d", i);
711 sprintf(numExc,
"N_EXC_%d", i);
712 sprintf(numIon,
"N_ION_%d", i);
713 sprintf(numPho,
"N_PHO_%d", i);
714 sprintf(numEle,
"N_ELE_%d", i);
715 NumExcitons = (G4int)aMaterialPropertiesTable->GetConstProperty(numExc);
716 NumIons = (G4int)aMaterialPropertiesTable->GetConstProperty(numIon);
717 sprintf(trackL,
"TRACK_%d", i);
718 sprintf(energy,
"ENRGY_%d", i);
719 sprintf(time00,
"TIME0_%d", i);
720 sprintf(time01,
"TIME1_%d", i);
721 delta = aMaterialPropertiesTable->GetConstProperty(trackL);
722 energ = aMaterialPropertiesTable->GetConstProperty(energy);
723 t0 = aMaterialPropertiesTable->GetConstProperty(time00);
724 t1 = aMaterialPropertiesTable->GetConstProperty(time01);
730 if ((delta < R0 && !fVeryHighEnergy) || z2 == z1 || fAlpha) {
731 if (z1 == 54 && ElectricField &&
732 Phase == kStateLiquid) {
734 abs(aParticle->GetPDGcode()) != 2112) {
735 ThomasImel = 0.056776 * pow(ElectricField, -0.11844);
737 ThomasImel = 0.057675 * pow(ElectricField, -0.49362);
742 ThomasImel = -0.15169 * pow((ElectricField + 215.25), 0.01811) + 0.20952;
745 if (ThomasImel > 0.19) ThomasImel = 0.19;
746 if (ThomasImel < 0.00) ThomasImel = 0.00;
748 if (Phase == kStateLiquid) {
749 if (z1 == 54) ThomasImel *= pow((Density / 2.84), 0.3);
750 if (z1 == 18) ThomasImel *= pow((Density /
Density_LAr), 0.3);
757 xi = (G4double(NumIons) / 4.) * ThomasImel;
758 if (InitialKinetEnergy == 9.4 * CLHEP::keV) {
759 G4double NumIonsEff =
760 1.066e7 * pow(t0 / CLHEP::ns, -1.303) * (0.17163 + 162.32 / (ElectricField + 191.39));
761 if (NumIonsEff > 1e6) NumIonsEff = 1e6;
762 xi = (G4double(NumIonsEff) / 4.) * ThomasImel;
764 if (fKr83m && ElectricField == 0) xi = (G4double(1.3 * NumIons) / 4.) * ThomasImel;
765 recombProb = 1 - log(1 + xi) / xi;
766 if (recombProb < 0) recombProb = 0;
767 if (recombProb > 1) recombProb = 1;
770 NumPhotons = NumExcitons +
BinomFluct(NumIons, recombProb);
771 NumElectrons = (NumExcitons + NumIons) - NumPhotons;
773 aMaterialPropertiesTable->AddConstProperty(numPho, NumPhotons);
774 aMaterialPropertiesTable->AddConstProperty(numEle, NumElectrons);
779 NumPhotons = (G4int)aMaterialPropertiesTable->GetConstProperty(numPho);
780 NumElectrons = (G4int)aMaterialPropertiesTable->GetConstProperty(numEle);
782 G4double FanoFactor = 0;
784 FanoFactor = 2575.9 * pow((ElectricField + 15.154), -0.64064) - 1.4707;
785 if (fKr83m) TotalEnergyDeposit = 4 * CLHEP::keV;
786 if ((dE / CLHEP::keV) <= 100 && ElectricField >= 0) {
787 G4double keVee = (TotalEnergyDeposit / (100. * CLHEP::keV));
789 FanoFactor *= -0.00075 + 0.4625 * keVee + 34.375 * pow(keVee, 2.);
791 FanoFactor *= 0.069554 + 1.7322 * keVee - .80215 * pow(keVee, 2.);
794 if (Phase == kStateGas && Density > 0.5)
795 FanoFactor = 0.42857 - 4.7857 * Density + 7.8571 * pow(Density, 2.);
796 if (FanoFactor <= 0 || fVeryHighEnergy) FanoFactor = 0;
797 NumQuanta = NumPhotons + NumElectrons;
798 if (z1 == 54 && FanoFactor)
800 G4int(floor(GaussGen.fire(NumElectrons, sqrt(FanoFactor * NumElectrons)) + 0.5));
801 NumPhotons = NumQuanta - NumElectrons;
802 if (NumElectrons <= 0) NumElectrons = 0;
810 NumElectrons = G4int(floor(NumElectrons *
phe_per_e + 0.5));
813 if (fKr83m || InitialKinetEnergy == 9.4 * CLHEP::keV) fKr83m += dE / CLHEP::keV;
814 if (fKr83m > 41) fKr83m = 0;
820 aMaterialPropertiesTable->AddConstProperty(numExc, 0);
821 aMaterialPropertiesTable->AddConstProperty(numIon, 0);
822 aMaterialPropertiesTable->AddConstProperty(numPho, 0);
823 aMaterialPropertiesTable->AddConstProperty(numEle, 0);
826 if (InitialKinetEnergy < MAX_ENE && InitialKinetEnergy >
MIN_ENE && !fMultipleScattering)
827 NumQuanta = NumPhotons + NumElectrons;
830 for (k = 0; k < NumQuanta; k++) {
831 G4double sampledEnergy;
832 std::unique_ptr<G4DynamicParticle> aQuantum;
835 G4double cost = 1. - 2. * UniformGen.fire();
836 G4double sint = std::sqrt((1. - cost) * (1. + cost));
837 G4double phi = CLHEP::twopi * UniformGen.fire();
838 G4double sinp = std::sin(phi);
839 G4double cosp = std::cos(phi);
840 G4double px = sint * cosp;
841 G4double py = sint * sinp;
845 G4ParticleMomentum photonMomentum(px, py, pz);
848 if (k < NumPhotons) {
850 G4double sx = cost * cosp;
851 G4double sy = cost * sinp;
853 G4ThreeVector photonPolarization(sx, sy, sz);
854 G4ThreeVector perp = photonMomentum.cross(photonPolarization);
855 phi = CLHEP::twopi * UniformGen.fire();
856 sinp = std::sin(phi);
857 cosp = std::cos(phi);
858 photonPolarization = cosp * photonPolarization + sinp * perp;
859 photonPolarization = photonPolarization.unit();
862 sampledEnergy = GaussGen.fire(PhotMean, PhotWidth);
864 std::make_unique<G4DynamicParticle>(G4OpticalPhoton::OpticalPhoton(), photonMomentum);
865 aQuantum->SetPolarization(
866 photonPolarization.x(), photonPolarization.y(), photonPolarization.z());
872 G4ParticleMomentum electronMomentum(0, 0, -FieldSign);
875 if (Phase == kStateGas) {
887 sampledEnergy = 1.38e-23 * (CLHEP::joule / CLHEP::kelvin) * Temperature;
892 aQuantum->SetKineticEnergy(sampledEnergy);
901 G4double aSecondaryTime = t0 + UniformGen.fire() * (t1 - t0) + evtStrt;
902 if (tau1 < 0) { tau1 = 0; }
903 if (tau3 < 0) { tau3 = 0; }
904 if (tauR < 0) { tauR = 0; }
905 if (aQuantum->GetDefinition()->GetParticleName() ==
"opticalphoton") {
906 if (
abs(z2 - z1) && !fAlpha &&
907 abs(aParticle->GetPDGcode()) != 2112) {
908 LET = (energ / CLHEP::MeV) / (delta / CLHEP::cm) * (1 / Density);
914 if (Phase == kStateLiquid && z1 == 54)
916 3.5 * ((1 + 0.41 * LET) / (0.18 * LET)) * CLHEP::ns * exp(-0.00900 * ElectricField);
923 SingTripRatioX = GaussGen.fire(0.17, 0.05);
924 SingTripRatioR = GaussGen.fire(0.8, 0.2);
927 0.2701 + 0.003379 * LET - 4.7338e-5 * pow(LET, 2.) + 8.1449e-6 * pow(LET, 3.);
928 SingTripRatioX = SingTripRatioR;
930 SingTripRatioX = GaussGen.fire(0.36, 0.06);
931 SingTripRatioR = GaussGen.fire(0.5, 0.2);
936 SingTripRatioR = GaussGen.fire(2.3, 0.51);
940 SingTripRatioR = (-0.065492 + 1.9996 * exp(-dE / CLHEP::MeV)) /
941 (1 + 0.082154 / pow(dE / CLHEP::MeV, 2.)) +
943 SingTripRatioX = SingTripRatioR;
948 SingTripRatioR = GaussGen.fire(7.8, 1.5);
949 if (z1 == 18) SingTripRatioR = 0.22218 * pow(energ / CLHEP::keV, 0.48211);
950 SingTripRatioX = SingTripRatioR;
955 if (k > NumExcitons) {
958 aSecondaryTime += tauR * (1. / UniformGen.fire() - 1);
959 if (UniformGen.fire() < SingTripRatioR / (1 + SingTripRatioR))
960 aSecondaryTime -= tau1 * log(UniformGen.fire());
962 aSecondaryTime -= tau3 * log(UniformGen.fire());
965 if (UniformGen.fire() < SingTripRatioX / (1 + SingTripRatioX))
966 aSecondaryTime -= tau1 * log(UniformGen.fire());
968 aSecondaryTime -= tau3 * log(UniformGen.fire());
972 G4double gainField = 12;
973 G4double tauTrap = 884.83 - 62.069 * gainField;
974 if (Phase == kStateLiquid) aSecondaryTime -= tauTrap * CLHEP::ns * log(UniformGen.fire());
983 std::string sx(xCoord);
984 std::string sy(yCoord);
985 std::string sz(zCoord);
986 x0[0] = aMaterialPropertiesTable->GetConstProperty(xCoord);
987 x0[1] = aMaterialPropertiesTable->GetConstProperty(yCoord);
988 x0[2] = aMaterialPropertiesTable->GetConstProperty(zCoord);
989 G4double
radius = sqrt(pow(x0[0], 2.) + pow(x0[1], 2.));
992 if (radius >=
R_TOL) {
993 if (x0[0] == 0) { x0[0] = 1 * CLHEP::nm; }
994 if (x0[1] == 0) { x0[1] = 1 * CLHEP::nm; }
996 phi = atan(x0[1] / x0[0]);
997 x0[0] = fabs(radius * cos(phi)) * ((fabs(x0[0])) / (x0[0]));
998 x0[1] = fabs(radius * sin(phi)) * ((fabs(x0[1])) / (x0[1]));
1001 G4ThreeVector aSecondaryPosition = x0;
1002 if (k >= NumPhotons &&
diffusion && ElectricField > 0) {
1003 G4double D_T = 64 * pow(1
e-3 * ElectricField, -.17);
1005 G4double D_L = 13.859 * pow(1
e-3 * ElectricField, -0.58559);
1007 if (Phase == kStateLiquid && z1 == 18) {
1008 D_T = 93.342 * pow(ElectricField / nDensity, 0.041322);
1011 if (Phase == kStateGas && z1 == 54) {
1012 D_L = 4.265 + 19097 / ElectricField - 1.7397e6 / pow(ElectricField, 2.) +
1013 1.2477e8 / pow(ElectricField, 3.);
1016 D_T *= CLHEP::cm2 / CLHEP::s;
1017 D_L *= CLHEP::cm2 / CLHEP::s;
1018 if (ElectricField < 100 && Phase == kStateLiquid) D_L = 8 * CLHEP::cm2 / CLHEP::s;
1019 G4double vDrift = sqrt((2 * sampledEnergy) / (
EMASS));
1020 if (
BORDER == 0) x0[2] = 0;
1021 G4double sigmaDT = sqrt(2 * D_T * fabs(
BORDER - x0[2]) / vDrift);
1022 G4double sigmaDL = sqrt(2 * D_L * fabs(
BORDER - x0[2]) / vDrift);
1023 G4double dr =
std::abs(GaussGen.fire(0., sigmaDT));
1024 phi = CLHEP::twopi * UniformGen.fire();
1025 aSecondaryPosition[0] += cos(phi) * dr;
1026 aSecondaryPosition[1] += sin(phi) * dr;
1027 aSecondaryPosition[2] += GaussGen.fire(0., sigmaDL);
1029 std::sqrt(std::pow(aSecondaryPosition[0], 2.) + std::pow(aSecondaryPosition[1], 2.));
1030 if (aSecondaryPosition[2] >=
BORDER && Phase == kStateLiquid) {
1033 if (aSecondaryPosition[2] <=
PMT && !GlobalFields) aSecondaryPosition[2] =
PMT +
R_TOL;
1037 if (aSecondaryTime < 0) aSecondaryTime = 0;
1048 aMaterialPropertiesTable->AddConstProperty(xCoord, 999 * CLHEP::km);
1049 aMaterialPropertiesTable->AddConstProperty(yCoord, 999 * CLHEP::km);
1050 aMaterialPropertiesTable->AddConstProperty(zCoord, 999 * CLHEP::km);
1051 aMaterialPropertiesTable->AddConstProperty(trackL, 0 * CLHEP::um);
1052 aMaterialPropertiesTable->AddConstProperty(energy, 0 * CLHEP::eV);
1053 aMaterialPropertiesTable->AddConstProperty(time00, DBL_MAX);
1054 aMaterialPropertiesTable->AddConstProperty(time01, -1 * CLHEP::ns);
1059 aMaterialPropertiesTable->AddConstProperty(
"TOTALNUM_INT_SITES", 0);
1060 aMaterialPropertiesTable->AddConstProperty(
"ENERGY_DEPOSIT_TOT", 0 * CLHEP::keV);
1061 aMaterialPropertiesTable->AddConstProperty(
"ENERGY_DEPOSIT_GOL", 0 * CLHEP::MeV);
1062 fExcitedNucleus =
false;
1078 std::cout <<
"WARNING: NestAlg::GetGasElectronDriftSpeed(G4double, G4double) " 1079 <<
"is not defined, returning bogus value of -999." << std::endl;
1086 G4double efieldinput,
1090 if (efieldinput < 0) efieldinput *= (-1);
1092 G4double onea = 144623.235704015, oneb = 850.812714257629, onec = 1192.87056676815,
1093 oned = -395969.575204061, onef = -355.484170008875, oneg = -227.266219627672,
1094 oneh = 223831.601257495, onei = 6.1778950907965, onej = 18.7831533426398,
1095 onek = -76132.6018884368;
1097 G4double twoa = 17486639.7118995, twob = -113.174284723134, twoc = 28.005913193763,
1098 twod = 167994210.094027, twof = -6766.42962575088, twog = 901.474643115395,
1099 twoh = -185240292.471665, twoi = -633.297790813084, twoj = 87.1756135457949;
1101 G4double thra = 10626463726.9833, thrb = 224025158.134792, thrc = 123254826.300172,
1102 thrd = -4563.5678061122, thrf = -1715.269592063, thrg = -694181.921834368,
1103 thrh = -50.9753281079838, thri = 58.3785811395493, thrj = 201512.080026704;
1104 G4double
y1 = 0,
y2 = 0,
f1 = 0,
f2 = 0,
f3 = 0, edrift = 0,
t1 = 0,
t2 = 0, slope = 0,
1108 f1 = onea / (1 + exp(-(efieldinput - oneb) / onec)) +
1109 oned / (1 + exp(-(efieldinput - onef) / oneg)) +
1110 oneh / (1 + exp(-(efieldinput - onei) / onej)) + onek;
1111 f2 = twoa / (1 + exp(-(efieldinput - twob) / twoc)) +
1112 twod / (1 + exp(-(efieldinput - twof) / twog)) +
1113 twoh / (1 + exp(-(efieldinput - twoi) / twoj));
1114 f3 = thra * exp(-thrb * efieldinput) + thrc * exp(-(pow(efieldinput - thrd, 2)) / (thrf * thrf)) +
1115 thrg * exp(-(pow(efieldinput - thrh, 2) / (thri * thri))) + thrj;
1117 if (efieldinput < 20 && efieldinput >= 0) {
1118 f1 = 2951 * efieldinput;
1119 f2 = 5312 * efieldinput;
1120 f3 = 7101 * efieldinput;
1123 if (tempinput < 200.0 && tempinput > 165.0) {
1129 if (tempinput < 230.0 && tempinput > 200.0) {
1135 if ((tempinput > 230.0 || tempinput < 165.0) && !Miller) {
1136 G4cout <<
"\nWARNING: TEMPERATURE OUT OF RANGE (165-230 K)\n";
1139 if (tempinput == 165.0)
1141 else if (tempinput == 200.0)
1143 else if (tempinput == 230.0)
1147 slope = (y1 -
y2) / (
t1 -
t2);
1148 intercept = y1 - slope *
t1;
1149 edrift = slope * tempinput + intercept;
1153 if (efieldinput <= 40.)
1154 edrift = -0.13274 + 0.041082 * efieldinput - 0.0006886 * pow(efieldinput, 2.) +
1155 5.5503e-6 * pow(efieldinput, 3.);
1157 edrift = 0.060774 * efieldinput / pow(1 + 0.11336 * pow(efieldinput, 0.5218), 2.);
1158 if (efieldinput >= 1e5) edrift = 2.7;
1159 if (efieldinput >= 100)
1160 edrift -= 0.017 * (tempinput - 163);
1162 edrift += 0.017 * (tempinput - 163);
1166 edrift = 1e5 * (.097384 * pow(log10(efieldinput), 3.0622) - .018614 * sqrt(efieldinput));
1167 if (edrift < 0) edrift = 0.;
1168 edrift = 0.5 *
EMASS * pow(edrift * CLHEP::cm / CLHEP::s, 2.);
1180 LET = 58.482 - 61.183 * log10(E) + 19.749 * pow(log10(E), 2) + 2.3101 * pow(log10(E), 3) -
1181 3.3469 * pow(log10(E), 4) + 0.96788 * pow(log10(E), 5) - 0.12619 * pow(log10(E), 6) +
1182 0.0065108 * pow(log10(E), 7);
1186 else if (E > 0 && E < 1)
1187 LET = 6.9463 + 815.98 * E - 4828 * pow(E, 2) + 17079 * pow(E, 3) - 36394 * pow(E, 4) +
1188 44553 * pow(E, 5) - 28659 * pow(E, 6) + 7483.8 * pow(E, 7);
1195 LET = 116.70 - 162.97 * log10(E) + 99.361 * pow(log10(E), 2) - 33.405 * pow(log10(E), 3) +
1196 6.5069 * pow(log10(E), 4) - 0.69334 * pow(log10(E), 5) + .031563 * pow(log10(E), 6);
1197 else if (E > 0 && E < 1)
1208 CLHEP::RandGauss GaussGen(
fEngine);
1209 CLHEP::RandFlat UniformGen(
fEngine);
1211 G4double
mean = N0 * prob;
1212 G4double sigma = sqrt(N0 * prob * (1 - prob));
1214 if (prob == 0.00)
return N1;
1215 if (prob == 1.00)
return N0;
1218 for (G4int i = 0; i < N0; i++) {
1219 if (UniformGen.fire() < prob) N1++;
1223 N1 = G4int(floor(GaussGen.fire(mean, sigma) + 0.5));
1225 if (N1 > N0) N1 = N0;
1246 for (G4int i = 0; i < 10000; i++) {
1247 sprintf(xCoord,
"POS_X_%d", i);
1248 sprintf(yCoord,
"POS_Y_%d", i);
1249 sprintf(zCoord,
"POS_Z_%d", i);
1250 nobleElementMat->AddConstProperty(xCoord, 999 * CLHEP::km);
1251 nobleElementMat->AddConstProperty(yCoord, 999 * CLHEP::km);
1252 nobleElementMat->AddConstProperty(zCoord, 999 * CLHEP::km);
1253 sprintf(numExc,
"N_EXC_%d", i);
1254 sprintf(numIon,
"N_ION_%d", i);
1255 sprintf(numPho,
"N_PHO_%d", i);
1256 sprintf(numEle,
"N_ELE_%d", i);
1257 nobleElementMat->AddConstProperty(numExc, 0);
1258 nobleElementMat->AddConstProperty(numIon, 0);
1259 nobleElementMat->AddConstProperty(numPho, 0);
1260 nobleElementMat->AddConstProperty(numEle, 0);
1261 sprintf(trackL,
"TRACK_%d", i);
1262 sprintf(energy,
"ENRGY_%d", i);
1263 sprintf(time00,
"TIME0_%d", i);
1264 sprintf(time01,
"TIME1_%d", i);
1265 nobleElementMat->AddConstProperty(trackL, 0 * CLHEP::um);
1266 nobleElementMat->AddConstProperty(energy, 0 * CLHEP::eV);
1267 nobleElementMat->AddConstProperty(time00, DBL_MAX);
1268 nobleElementMat->AddConstProperty(time01, -1 * CLHEP::ns);
1274 nobleElementMat->AddConstProperty(
"TOTALNUM_INT_SITES", 0);
1275 nobleElementMat->AddConstProperty(
"ENERGY_DEPOSIT_TOT", 0 * CLHEP::keV);
1276 nobleElementMat->AddConstProperty(
"ENERGY_DEPOSIT_GOL", 0 * CLHEP::MeV);
1286 G4double a_0 = 5.29e-11 * CLHEP::m;
1287 G4double a = 0.626 * a_0 * pow(Z, (-1. / 3.));
1288 G4double epsilon_0 = 8.854e-12 * (CLHEP::farad / CLHEP::m);
1290 (a * E * 2. * CLHEP::twopi * epsilon_0) / (2 * pow(CLHEP::eplus, 2.) * pow(Z, 2.));
1291 G4double zeta_0 = pow(Z, (1. / 6.));
1292 G4double m_N = A * 1.66e-27 * CLHEP::kg;
1293 G4double hbar = 6.582e-16 * CLHEP::eV * CLHEP::s;
1297 G4double s_n = log(1 + 1.1383 * epsilon) /
1298 (2. * (epsilon + 0.01321 * pow(epsilon, 0.21226) + 0.19593 * sqrt(epsilon)));
1300 (a_0 * zeta_0 / a) * hbar *
1301 sqrt(8 * epsilon * 2. * CLHEP::twopi * epsilon_0 / (a * m_N * pow(CLHEP::eplus, 2.)));
1302 return 1.38e5 * 0.5 * (1 + tanh(50 * epsilon - 0.25)) * epsilon * (s_e / s_n);
code to link reconstructed objects back to the MC truth information
Float_t y1[n_points_granero]
G4double GetLiquidElectronDriftSpeed(double T, double F, G4bool M, G4int Z)
Float_t x1[n_points_granero]
NestAlg(CLHEP::HepRandomEngine &engine)
constexpr auto abs(T v)
Returns the absolute value of the argument.
double fYieldFactor
turns scint. on/off
G4VParticleChange fParticleChange
pointer to G4VParticleChange
int fNumIonElectrons
number of ionization electrons produced by step
kilovolt_as<> kilovolt
Type of potential stored in kilovolt, in double precision.
G4int BinomFluct(G4int N0, G4double prob)
Float_t y2[n_points_geant4]
int fNumScintPhotons
number of photons produced by the step
auto counter(T begin, T end)
Returns an object to iterate values from begin to end in a range-for loop.
static G4ThermalElectron * ThermalElectron()
const G4VParticleChange & CalculateIonizationAndScintillation(G4Track const &aTrack, G4Step const &aStep)
G4double UnivScreenFunc(G4double E, G4double Z, G4double A)
G4double GetGasElectronDriftSpeed(G4double efieldinput, G4double density)
CLHEP::HepRandomEngine & fEngine
random engine
G4double CalculateElectronLET(G4double E, G4int Z)
double mean(const std::vector< short > &wf, size_t start, size_t nsample)
void InitMatPropValues(G4MaterialPropertiesTable *nobleElementMat, int z)
std::map< int, bool > fElementPropInit
double fEnergyDep
energy deposited by the step