| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 ///////////////////////////////////////////////////////////////////////
28 // UCN BoundaryProcess Class Implementation
29 ///////////////////////////////////////////////////////////////////////
30 //
31 // File: G4UCNBoundaryProcess.cc
32 // Description: Discrete Process -- Boundary Process of UCN
33 // Version: 1.0
34 // Created: 2014-05-12
35 // Author: Peter Gumplinger
36 // adopted from Geant4UCN by Peter Fierlinger (9.9.04) and
37 // Marcin Kuzniak (21.4.06)
38 // 1/v energy dependent absorption cross section
39 // inside materials
40 // Updated: 2007 Extensions for the microroughness model by Stefan Heule
41 //
42 // mail: gum@triumf.ca
43 //
44 ///////////////////////////////////////////////////////////////////////
45
46 #include "G4UCNProcessSubType.hh"
47
48
49 #include "G4UCNBoundaryProcess.hh"
50 #include "G4UCNBoundaryProcessMessenger.hh"
51
52 #include "G4GeometryTolerance.hh"
53
54 #include "G4StepPoint.hh"
55 #include "G4ParticleDefinition.hh"
56
57 #include "G4UCNMaterialPropertiesTable.hh"
58
59 #include "G4TransportationManager.hh"
60 #include "G4ParallelWorldProcess.hh"
61
62 #include "G4VSensitiveDetector.hh"
63
64 #include "G4SystemOfUnits.hh"
65 #include "G4PhysicalConstants.hh"
66
67 G4UCNBoundaryProcess::G4UCNBoundaryProcess(const G4String& processName,
68 G4ProcessType type)
69 : G4VDiscreteProcess(processName, type)
70 {
71
72 if (verboseLevel > 0) G4cout << GetProcessName() << " is created " << G4endl;
73
74 SetProcessSubType(fUCNBoundary);
75
76 theStatus = Undefined;
77
78 fMessenger = new G4UCNBoundaryProcessMessenger(this);
79
80 neV = 1.0e-9*eV;
81
82 kCarTolerance = G4GeometryTolerance::GetInstance()
83 ->GetSurfaceTolerance();
84
85 Material1 = NULL;
86 Material2 = NULL;
87
88 aMaterialPropertiesTable1 = NULL;
89 aMaterialPropertiesTable2 = NULL;
90
91 UseMicroRoughnessReflection = false;
92 DoMicroRoughnessReflection = false;
93
94 nNoMPT = nNoMRT = nNoMRCondition = 0;
95 nAbsorption = nEzero = nFlip = 0;
96 aSpecularReflection = bSpecularReflection = 0;
97 bLambertianReflection = 0;
98 aMRDiffuseReflection = bMRDiffuseReflection = 0;
99 nSnellTransmit = mSnellTransmit = 0;
100 aMRDiffuseTransmit = 0;
101 ftheta_o = fphi_o = 0;
102 }
103
104 G4UCNBoundaryProcess::~G4UCNBoundaryProcess()
105 {
106 delete fMessenger;
107 }
108
109 G4VParticleChange*
110 G4UCNBoundaryProcess::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
111 {
112 aParticleChange.Initialize(aTrack);
113 aParticleChange.ProposeVelocity(aTrack.GetVelocity());
114
115 // Get hyperStep from G4ParallelWorldProcess
116 // NOTE: PostSetpDoIt of this process should be
117 // invoked after G4ParallelWorldProcess!
118
119 const G4Step* pStep = &aStep;
120
121 const G4Step* hStep = G4ParallelWorldProcess::GetHyperStep();
122
123 if (hStep) pStep = hStep;
124
125 G4bool isOnBoundary =
126 (pStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary);
127
128 if (isOnBoundary) {
129 Material1 = pStep->GetPreStepPoint()->GetMaterial();
130 Material2 = pStep->GetPostStepPoint()->GetMaterial();
131 } else {
132 theStatus = NotAtBoundary;
133 if ( verboseLevel > 1 ) BoundaryProcessVerbose();
134 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
135 }
136
137 if (aTrack.GetStepLength()<=kCarTolerance/2) {
138 theStatus = StepTooSmall;
139 if ( verboseLevel > 0 ) BoundaryProcessVerbose();
140 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
141 }
142
143 aMaterialPropertiesTable1 = (G4UCNMaterialPropertiesTable*)Material1->
144 GetMaterialPropertiesTable();
145 aMaterialPropertiesTable2 = (G4UCNMaterialPropertiesTable*)Material2->
146 GetMaterialPropertiesTable();
147
148 G4String volnam1 = pStep->GetPreStepPoint() ->GetPhysicalVolume()->GetName();
149 G4String volnam2 = pStep->GetPostStepPoint()->GetPhysicalVolume()->GetName();
150
151 if (verboseLevel > 0) {
152 G4cout << " UCN at Boundary! " << G4endl;
153 G4cout << " vol1: " << volnam1 << ", vol2: " << volnam2 << G4endl;
154 G4cout << " Ekin: " << aTrack.GetKineticEnergy()/neV <<"neV"<< G4endl;
155 G4cout << " MomDir: " << aTrack.GetMomentumDirection() << G4endl;
156 }
157
158 if (Material1 == Material2) {
159 theStatus = SameMaterial;
160 if ( verboseLevel > 0 ) BoundaryProcessVerbose();
161 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
162 }
163
164 G4ThreeVector theGlobalPoint = pStep->GetPostStepPoint()->GetPosition();
165
166 G4bool valid;
167 // Use the new method for Exit Normal in global coordinates,
168 // which provides the normal more reliably.
169
170 // ID of Navigator which limits step
171
172 G4int hNavId = G4ParallelWorldProcess::GetHypNavigatorID();
173 std::vector<G4Navigator*>::iterator iNav =
174 G4TransportationManager::GetTransportationManager()->
175 GetActiveNavigatorsIterator();
176
177 G4ThreeVector theGlobalNormal =
178 (iNav[hNavId])->GetGlobalExitNormal(theGlobalPoint,&valid);
179
180 if (valid) {
181 theGlobalNormal = -theGlobalNormal;
182 }
183 else
184 {
185 G4ExceptionDescription ed;
186 ed << " G4UCNBoundaryProcess/PostStepDoIt(): "
187 << " The Navigator reports that it returned an invalid normal"
188 << G4endl;
189 G4Exception("G4UCNBoundaryProcess::PostStepDoIt", "UCNBoun01",
190 EventMustBeAborted,ed,
191 "Invalid Surface Normal - Geometry must return valid surface normal");
192 }
193
194 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
195
196 G4ThreeVector OldMomentum = aParticle->GetMomentumDirection();
197
198 if (OldMomentum * theGlobalNormal > 0.0) {
199 #ifdef G4OPTICAL_DEBUG
200 G4ExceptionDescription ed;
201 ed << " G4UCNBoundaryProcess/PostStepDoIt(): "
202 << " theGlobalNormal points in a wrong direction. "
203 << G4endl;
204 ed << " The momentum of the photon arriving at interface (oldMomentum)"
205 << " must exit the volume cross in the step. " << G4endl;
206 ed << " So it MUST have dot < 0 with the normal that Exits the new volume (globalNormal)." << G4endl;
207 ed << " >> The dot product of oldMomentum and global Normal is " << OldMomentum*theGlobalNormal << G4endl;
208 ed << " Old Momentum (during step) = " << OldMomentum << G4endl;
209 ed << " Global Normal (Exiting New Vol) = " << theGlobalNormal << G4endl;
210 ed << G4endl;
211 G4Exception("G4UCNBoundaryProcess::PostStepDoIt", "UCNBoun02",
212 EventMustBeAborted, // Or JustWarning to see if it happens repeatedbly on one ray
213 ed,
214 "Invalid Surface Normal - Geometry must return valid surface normal pointing in the right direction");
215 #else
216 theGlobalNormal = -theGlobalNormal;
217 #endif
218 }
219
220 G4ThreeVector theNeutronMomentum = aTrack.GetMomentum();
221
222 G4double theMomentumNormal = theNeutronMomentum*theGlobalNormal;
223 G4double theVelocityNormal = aTrack.GetVelocity() *
224 (OldMomentum * theGlobalNormal);
225
226 G4double Enormal = theMomentumNormal*theMomentumNormal/2./neutron_mass_c2;
227 G4double Energy = aTrack.GetKineticEnergy();
228
229 G4double FermiPot2 = 0.;
230 G4double pDiffuse = 0.;
231 G4double pSpinFlip = 0.;
232 G4double pUpScatter = 0.;
233
234 if (aMaterialPropertiesTable2) {
235 FermiPot2 = aMaterialPropertiesTable2->GetConstProperty("FERMIPOT")*neV;
236 pDiffuse = aMaterialPropertiesTable2->GetConstProperty("DIFFUSION");
237 pSpinFlip = aMaterialPropertiesTable2->GetConstProperty("SPINFLIP");
238 pUpScatter = aMaterialPropertiesTable2->GetConstProperty("LOSS");
239 }
240
241 G4double FermiPot1 = 0.;
242 if (aMaterialPropertiesTable1)
243 FermiPot1 = aMaterialPropertiesTable1->GetConstProperty("FERMIPOT")*neV;
244
245 G4double FermiPotDiff = FermiPot2 - FermiPot1;
246
247 if ( verboseLevel > 1 )
248 G4cout << "UCNBoundaryProcess: new FermiPot: " << FermiPot2/neV
249 << "neV, old FermiPot:" << FermiPot1/neV << "neV" << G4endl;
250
251 // Use microroughness diffuse reflection behavior if activated
252
253 DoMicroRoughnessReflection = UseMicroRoughnessReflection;
254
255 G4double theta_i = 0.;
256
257 if (!aMaterialPropertiesTable2) {
258
259 nNoMPT++;
260 theStatus = NoMPT;
261 if ( verboseLevel > 0 ) BoundaryProcessVerbose();
262 DoMicroRoughnessReflection = false;
263
264 } else {
265
266 if (!aMaterialPropertiesTable2->GetMicroRoughnessTable()) {
267
268 nNoMRT++;
269 theStatus = NoMRT;
270 if ( verboseLevel > 0 ) BoundaryProcessVerbose();
271
272 DoMicroRoughnessReflection = false;
273 }
274
275 // Angle theta_in between surface and momentum direction,
276 // Phi_in is defined to be 0
277
278 theta_i = OldMomentum.angle(-theGlobalNormal);
279
280 // Checks the MR-conditions
281
282 if (!aMaterialPropertiesTable2->
283 ConditionsValid(Energy, FermiPotDiff, theta_i)) {
284
285 nNoMRCondition++;
286 theStatus = NoMRCondition;
287 if ( verboseLevel > 0 ) BoundaryProcessVerbose();
288
289 DoMicroRoughnessReflection = false;
290 }
291 }
292
293 G4double MRpDiffuse = 0.;
294 G4double MRpDiffuseTrans = 0.;
295
296 // If microroughness is available and active for material in question
297
298 if (DoMicroRoughnessReflection) {
299
300 // Integral probability for non-specular reflection with microroughness
301
302 MRpDiffuse = aMaterialPropertiesTable2->
303 GetMRIntProbability(theta_i, Energy);
304
305 // Integral probability for non-specular transmission with microroughness
306
307 MRpDiffuseTrans = aMaterialPropertiesTable2->
308 GetMRIntTransProbability(theta_i, Energy);
309
310 if ( verboseLevel > 1 ) {
311 G4cout << "theta: " << theta_i/degree << "degree" << G4endl;
312 G4cout << "Energy: " << Energy/neV << "neV"
313 << ", Enormal: " << Enormal/neV << "neV" << G4endl;
314 G4cout << "FermiPotDiff: " << FermiPotDiff/neV << "neV " << G4endl;
315 G4cout << "Reflection_prob: " << MRpDiffuse
316 << ", Transmission_prob: " << MRpDiffuseTrans << G4endl;
317 }
318 }
319
320 if (!High(Enormal, FermiPotDiff)) {
321
322 // Below critical velocity
323
324 if (verboseLevel > 0) G4cout << "G4UCNBoundaryProcess -> BELOW critical velocity" << G4endl;
325
326 // Loss on reflection
327
328 if (Loss(pUpScatter, theVelocityNormal, FermiPotDiff)) {
329
330 // kill it.
331 aParticleChange.ProposeTrackStatus(fStopAndKill);
332 aParticleChange.ProposeLocalEnergyDeposit(Energy);
333
334 nAbsorption++;
335 theStatus = Absorption;
336 if ( verboseLevel > 0 ) BoundaryProcessVerbose();
337
338 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
339 }
340
341 // spinflips
342
343 if (SpinFlip(pSpinFlip)) {
344 nFlip++;
345 theStatus = Flip;
346 if ( verboseLevel > 0 ) BoundaryProcessVerbose();
347
348 G4ThreeVector NewPolarization = -1. * aParticle->GetPolarization();
349 aParticleChange.ProposePolarization(NewPolarization);
350 }
351
352 // Reflect from surface
353
354 G4ThreeVector NewMomentum;
355
356 // If microroughness is available and active - do non-specular reflection
357
358 if (DoMicroRoughnessReflection)
359 NewMomentum = MRreflect(MRpDiffuse, OldMomentum, theGlobalNormal,
360 Energy, FermiPotDiff);
361 else
362
363 // Else do it with the Lambert model as implemented by Peter Fierlinger
364
365 NewMomentum = Reflect(pDiffuse, OldMomentum, theGlobalNormal);
366
367 aParticleChange.ProposeMomentumDirection(NewMomentum);
368
369 } else {
370
371 // Above critical velocity
372
373 if (verboseLevel > 0) G4cout << "G4UCNBoundaryProcess -> ABOVE critical velocity" << G4endl;
374
375 // If it is faster than the criticial velocity,
376 // there is a probability to be still reflected.
377 // This formula is (only) valid for low loss materials
378
379 // If microroughness available and active, do reflection with it
380
381 G4ThreeVector NewMomentum;
382
383 if (DoMicroRoughnessReflection) {
384
385 G4double Enew;
386
387 NewMomentum =
388 MRreflectHigh(MRpDiffuse, MRpDiffuseTrans, 0., OldMomentum,
389 theGlobalNormal, Energy, FermiPotDiff, Enew);
390
391 if (Enew == 0.) {
392 aParticleChange.ProposeTrackStatus(fStopAndKill);
393 aParticleChange.ProposeLocalEnergyDeposit(Energy);
394 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
395 } else {
396 aParticleChange.ProposeEnergy(Enew);
397 aParticleChange.ProposeMomentumDirection(NewMomentum);
398 aParticleChange.ProposeVelocity(std::sqrt(2*Enew/neutron_mass_c2)*c_light);
399 aParticleChange.ProposeLocalEnergyDeposit(Energy-Enew);
400 }
401
402 } else {
403
404 G4double reflectivity = Reflectivity(FermiPotDiff, Enormal);
405
406 if ( verboseLevel > 1 ) G4cout << "UCNBoundaryProcess: reflectivity "
407 << reflectivity << G4endl;
408
409 if (G4UniformRand() < reflectivity) {
410
411 // Reflect from surface
412
413 NewMomentum = Reflect(pDiffuse, OldMomentum, theGlobalNormal);
414 aParticleChange.ProposeMomentumDirection(NewMomentum);
415
416 } else {
417
418 // --- Transmission because it is faster than the critical velocity
419
420 G4double Enew = Transmit(FermiPotDiff, Energy);
421
422 // --- Change of the normal momentum component
423 // p = sqrt(2*m*Ekin)
424
425 G4double mass = -std::sqrt(theMomentumNormal*theMomentumNormal -
426 neutron_mass_c2*2.*FermiPotDiff);
427
428 // --- Momentum direction in new media
429
430 NewMomentum =
431 theNeutronMomentum - (theMomentumNormal-mass)*theGlobalNormal;
432
433 nSnellTransmit++;
434 theStatus = SnellTransmit;
435 if ( verboseLevel > 0 ) BoundaryProcessVerbose();
436
437 aParticleChange.ProposeEnergy(Enew);
438 aParticleChange.ProposeMomentumDirection(NewMomentum.unit());
439 aParticleChange.ProposeVelocity(std::sqrt(2*Enew/neutron_mass_c2)*c_light);
440 aParticleChange.ProposeLocalEnergyDeposit(Energy-Enew);
441
442 if (verboseLevel > 1) {
443 G4cout << "Energy: " << Energy/neV << "neV, Enormal: "
444 << Enormal/neV << "neV, fpdiff " << FermiPotDiff/neV
445 << "neV, Enew " << Enew/neV << "neV" << G4endl;
446 G4cout << "UCNBoundaryProcess: Transmit and set the new Energy "
447 << aParticleChange.GetEnergy()/neV << "neV" << G4endl;
448 }
449 }
450 }
451 }
452
453 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
454 }
455
456 G4double G4UCNBoundaryProcess::GetMeanFreePath(const G4Track&,
457 G4double ,
458 G4ForceCondition* condition)
459 {
460 *condition = Forced;
461
462 return DBL_MAX;
463 }
464
465 G4bool G4UCNBoundaryProcess::Loss(G4double pUpScatter,
466 G4double theVelocityNormal,
467 G4double FermiPot)
468 {
469 // The surface roughness is not taken into account here.
470 // One could use e.g. ultracold neutrons, R. Golub, p.35,
471 // where mu is increased by roughness parameters sigma and
472 // omega, which are related to the height and the distance of
473 // "disturbances" on the surface
474
475 G4double vBound = std::sqrt(2.*FermiPot/neutron_mass_c2*c_squared);
476 G4double vRatio = theVelocityNormal/vBound;
477
478 G4double pLoss = (2*pUpScatter*vRatio)/(std::sqrt(1-(vRatio*vRatio)));
479
480 // Check, if enhancement for surface roughness should be done
481
482 if (DoMicroRoughnessReflection) {
483 if (aMaterialPropertiesTable2) {
484 const G4double hdm = hbar_Planck*c_squared/neutron_mass_c2;
485 G4double b = aMaterialPropertiesTable2->GetRMS();
486 G4double w = aMaterialPropertiesTable2->GetCorrLen();
487
488 // cf. Golub's book p. 35, eq. 2.103
489
490 pLoss *= std::sqrt(1+2*b*b*vBound*vBound/
491 (hdm*hdm+0.85*hdm*vBound*w+2*vBound*vBound*w*w));
492 }
493 }
494
495 return (G4UniformRand() <= std::fabs(pLoss));
496 }
497
498 G4bool G4UCNBoundaryProcess::SpinFlip(G4double pSpinFlip)
499 {
500 return (G4UniformRand() <= pSpinFlip);
501 }
502
503 G4double G4UCNBoundaryProcess::Reflectivity(G4double FermiPot, G4double Enormal)
504 {
505 G4double r = (std::sqrt(Enormal) - std::sqrt(Enormal - FermiPot)) /
506 (std::sqrt(Enormal) + std::sqrt(Enormal - FermiPot));
507
508 return r*r;
509 }
510
511 G4ThreeVector G4UCNBoundaryProcess::Reflect(G4double pDiffuse,
512 G4ThreeVector OldMomentum,
513 G4ThreeVector Normal)
514 {
515 G4double PdotN = OldMomentum * Normal;
516
517 G4ThreeVector NewMomentum = OldMomentum - (2.*PdotN)*Normal;
518 NewMomentum.unit();
519
520 // Reflect diffuse
521
522 if (NewMomentum == OldMomentum || G4UniformRand() < pDiffuse) {
523
524 NewMomentum = LDiffRefl(Normal);
525
526 bLambertianReflection++;
527 theStatus = LambertianReflection;
528 if ( verboseLevel > 0 ) BoundaryProcessVerbose();
529
530 return NewMomentum;
531 }
532
533 // Reflect specular
534
535 bSpecularReflection++;
536 theStatus = SpecularReflection;
537 if ( verboseLevel > 0 ) BoundaryProcessVerbose();
538
539 return NewMomentum;
540 }
541
542 G4ThreeVector G4UCNBoundaryProcess::MRreflect(G4double pDiffuse,
543 G4ThreeVector OldMomentum,
544 G4ThreeVector Normal,
545 G4double Energy,
546 G4double FermiPot)
547 {
548 // Only for Enormal <= VFermi
549
550 G4ThreeVector NewMomentum;
551
552 // Checks if the reflection should be non-specular
553
554 if (G4UniformRand() <= pDiffuse) {
555
556 // Reflect diffuse
557
558 // Determines the angles of the non-specular reflection
559
560 NewMomentum =
561 MRDiffRefl(Normal, Energy, FermiPot, OldMomentum, pDiffuse);
562
563 bMRDiffuseReflection++;
564 theStatus = MRDiffuseReflection;
565 if ( verboseLevel > 0 ) BoundaryProcessVerbose();
566
567 return NewMomentum;
568
569 } else {
570
571 // Reflect specular
572
573 G4double PdotN = OldMomentum * Normal;
574
575 NewMomentum = OldMomentum - (2.*PdotN)*Normal;
576 NewMomentum.unit();
577
578 bSpecularReflection++;
579 theStatus = SpecularReflection;
580 if ( verboseLevel > 0 ) BoundaryProcessVerbose();
581
582 return NewMomentum;
583 }
584 }
585
586 G4ThreeVector G4UCNBoundaryProcess::MRreflectHigh(G4double pDiffuse,
587 G4double pDiffuseTrans,
588 G4double pLoss,
589 G4ThreeVector OldMomentum,
590 G4ThreeVector Normal,
591 G4double Energy,
592 G4double FermiPot,
593 G4double &Enew)
594 {
595 // Only for Enormal > VFermi
596
597 G4double costheta = OldMomentum*Normal;
598
599 G4double Enormal = Energy * (costheta*costheta);
600
601 // G4double pSpecular = Reflectivity(Enormal,FermiPot)*
602 G4double pSpecular = Reflectivity(FermiPot,Enormal)*
603 (1.-pDiffuse-pDiffuseTrans-pLoss);
604
605 G4ThreeVector NewMomentum;
606
607 G4double decide = G4UniformRand();
608
609 if (decide < pSpecular) {
610
611 // Reflect specularly
612
613 G4double PdotN = OldMomentum * Normal;
614 NewMomentum = OldMomentum - (2.*PdotN)*Normal;
615 NewMomentum.unit();
616
617 Enew = Energy;
618
619 aSpecularReflection++;
620 theStatus = SpecularReflection;
621 if ( verboseLevel ) BoundaryProcessVerbose();
622
623 return NewMomentum;
624 }
625
626 if (decide < pSpecular + pDiffuse) {
627
628 // Reflect diffusely
629
630 // Determines the angles of the non-specular reflection
631
632 NewMomentum =
633 MRDiffRefl(Normal, Energy, FermiPot, OldMomentum, pDiffuse);
634
635 if (verboseLevel > 0) G4cout << "Diffuse normal " << Normal
636 << ", " << NewMomentum << G4endl;
637 Enew = Energy;
638
639 aMRDiffuseReflection++;
640 theStatus = MRDiffuseReflection;
641 if ( verboseLevel ) BoundaryProcessVerbose();
642
643 return NewMomentum;
644 }
645
646 if (decide < pSpecular + pDiffuse + pDiffuseTrans) {
647
648 // Transmit diffusely
649
650 // Determines the angles of the non-specular transmission
651
652 NewMomentum =
653 MRDiffTrans(Normal, Energy, FermiPot, OldMomentum, pDiffuseTrans);
654
655 Enew = Energy - FermiPot;
656
657 aMRDiffuseTransmit++;
658 theStatus = MRDiffuseTransmit;
659 if ( verboseLevel ) BoundaryProcessVerbose();
660
661 return NewMomentum;
662 }
663
664 if (decide < pSpecular + pDiffuse + pDiffuseTrans + pLoss) {
665
666 // Loss
667
668 Enew = 0.;
669
670 nEzero++;
671 theStatus = Ezero;
672 if ( verboseLevel > 0 ) BoundaryProcessVerbose();
673
674 return NewMomentum;
675 }
676
677 // last case: Refractive transmission
678
679 Enew = Energy - FermiPot;
680
681 G4double palt = std::sqrt(2*neutron_mass_c2/c_squared*Energy);
682 G4double produ = OldMomentum * Normal;
683
684 NewMomentum = palt*OldMomentum-
685 (palt*produ+std::sqrt(palt*palt*produ*produ-2*neutron_mass_c2/
686 c_squared*FermiPot))*Normal;
687
688 mSnellTransmit++;
689 theStatus = SnellTransmit;
690 if ( verboseLevel > 0 ) BoundaryProcessVerbose();
691
692 return NewMomentum.unit();
693 }
694
695 G4ThreeVector G4UCNBoundaryProcess::MRDiffRefl(G4ThreeVector Normal,
696 G4double Energy,
697 G4double FermiPot,
698 G4ThreeVector OldMomentum,
699 G4double pDiffuse)
700 {
701 G4bool accepted = false;
702
703 G4double theta_o, phi_o;
704
705 // Polar angle of incidence
706
707 G4double theta_i = OldMomentum.polarAngle(-Normal);
708
709 // Azimuthal angle of incidence
710
711 // G4double phi_i = -OldMomentum.azimAngle(-Normal);
712
713 // accept-reject method for MR-reflection
714
715 G4int count = 0;
716 while (!accepted) {
717 theta_o = G4UniformRand()*pi/2;
718 phi_o = G4UniformRand()*pi*2-pi;
719 // Box over distribution is increased by 50% to ensure no value is above
720 if (1.5*G4UniformRand()*
721 aMaterialPropertiesTable2->
722 GetMRMaxProbability(theta_i, Energy)/pDiffuse <=
723 aMaterialPropertiesTable2->
724 GetMRProbability(theta_i,Energy,FermiPot,theta_o,phi_o)/pDiffuse)
725
726 accepted = true;
727
728 // For the case that the box is nevertheless exceeded
729
730 if (aMaterialPropertiesTable2->
731 GetMRProbability(theta_i, Energy, FermiPot, theta_o, phi_o)/
732 (1.5*aMaterialPropertiesTable2->
733 GetMRMaxProbability(theta_i, Energy)) > 1) {
734 G4cout << "MRMax Wahrscheinlichkeitsueberschreitung!" << G4endl;
735 G4cout << aMaterialPropertiesTable2->
736 GetMRProbability(theta_i, Energy, FermiPot, theta_o, phi_o)/
737 (1.5*aMaterialPropertiesTable2->
738 GetMRMaxProbability(theta_i, Energy)) << G4endl;
739 aMaterialPropertiesTable2->
740 SetMRMaxProbability(theta_i, Energy,
741 aMaterialPropertiesTable2->
742 GetMRProbability(theta_i, Energy,
743 FermiPot, theta_o, phi_o));
744 }
745 // Loop checking, 31-Aug-2015, Vladimir Ivanchenko
746 if(++count > 10000) { accepted = true; }
747 }
748
749 // Creates vector in the local coordinate system of the reflection
750
751 G4ThreeVector localmomentum;
752 localmomentum.setRThetaPhi(1., theta_o, phi_o);
753
754 ftheta_o = theta_o;
755 fphi_o = phi_o;
756
757 // Get coordinate transform matrix
758
759 G4RotationMatrix TransCoord =
760 GetCoordinateTransformMatrix(Normal, OldMomentum);
761
762 //transfer to the coordinates of the global system
763
764 G4ThreeVector momentum = TransCoord*localmomentum;
765
766 //momentum.rotateUz(Normal);
767
768 if (momentum * Normal<0) {
769 momentum*=-1;
770 // something has gone wrong...
771 G4cout << "G4UCNBoundaryProcess::MRDiffRefl: !" << G4endl;
772 }
773
774 return momentum.unit();
775 }
776
777 G4ThreeVector G4UCNBoundaryProcess::MRDiffTrans(G4ThreeVector Normal,
778 G4double Energy,
779 G4double FermiPot,
780 G4ThreeVector OldMomentum,
781 G4double pDiffuseTrans)
782 {
783 G4bool accepted = false;
784
785 G4double theta_o, phi_o;
786
787 // Polar angle of incidence
788
789 G4double theta_i = OldMomentum.polarAngle(-Normal);
790
791 // azimuthal angle of incidence
792
793 // G4double phi_i = -OldMomentum.azimAngle(-Normal);
794
795 G4int count = 0;
796 while (!accepted) {
797 theta_o = G4UniformRand()*pi/2;
798 phi_o = G4UniformRand()*pi*2-pi;
799
800 // Box over distribution is increased by 50% to ensure no value is above
801
802 if (1.5*G4UniformRand()*
803 aMaterialPropertiesTable2->
804 GetMRMaxTransProbability(theta_i, Energy)/pDiffuseTrans <=
805 aMaterialPropertiesTable2->
806 GetMRTransProbability(theta_i,Energy,FermiPot,theta_o,phi_o)/
807 pDiffuseTrans)
808
809 accepted=true;
810
811 // For the case that the box is nevertheless exceeded
812
813 if(aMaterialPropertiesTable2->
814 GetMRTransProbability(theta_i, Energy, FermiPot, theta_o, phi_o)/
815 (1.5*aMaterialPropertiesTable2->
816 GetMRMaxTransProbability(theta_i, Energy)) > 1) {
817 G4cout << "MRMaxTrans Wahrscheinlichkeitsueberschreitung!" << G4endl;
818 G4cout << aMaterialPropertiesTable2->
819 GetMRTransProbability(theta_i, Energy, FermiPot, theta_o, phi_o)/
820 (1.5*aMaterialPropertiesTable2->
821 GetMRMaxTransProbability(theta_i, Energy)) << G4endl;
822 aMaterialPropertiesTable2->
823 SetMRMaxTransProbability(theta_i, Energy,
824 aMaterialPropertiesTable2->
825 GetMRTransProbability(theta_i, Energy,
826 FermiPot, theta_o, phi_o));
827 }
828 // Loop checking, 31-Aug-2015, Vladimir Ivanchenko
829 if(++count > 10000) { accepted = true; }
830 }
831
832 // Creates vector in the local coordinate system of the reflection
833
834 G4ThreeVector localmomentum;
835 localmomentum.setRThetaPhi(1., pi-theta_o, phi_o);
836
837 // Get coordinate transform matrix
838
839 G4RotationMatrix TransCoord =
840 GetCoordinateTransformMatrix(Normal, OldMomentum);
841
842 // Transfer to the coordinates of the global system
843
844 G4ThreeVector momentum = TransCoord*localmomentum;
845
846 if (momentum*Normal<0) {
847 // something has gone wrong...
848 momentum*=-1;
849 G4cout << "G4UCNBoundaryProcess::MRDiffTrans: !" << G4endl;
850 }
851
852 return momentum.unit();
853 }
854
855 G4double G4UCNBoundaryProcess::Transmit(G4double FermiPot, G4double Energy)
856 {
857 return (Energy - FermiPot);
858 }
859
860 G4ThreeVector G4UCNBoundaryProcess::LDiffRefl(G4ThreeVector Normal)
861 {
862 G4ThreeVector momentum;
863
864 // cosine distribution - Lambert's law
865
866 momentum.setRThetaPhi(1., std::acos(std::sqrt(G4UniformRand())), 2.*pi*G4UniformRand());
867 momentum.rotateUz(Normal);
868
869 if (momentum*Normal < 0) {
870 momentum*=-1;
871 G4cout << "G4UCNBoundaryProcess::LDiffRefl: !" << G4endl;
872 }
873
874 return momentum.unit();
875 }
876
877 G4RotationMatrix G4UCNBoundaryProcess::
878 GetCoordinateTransformMatrix(G4ThreeVector Normal,
879 G4ThreeVector direction)
880 {
881 // Definition of the local coordinate system (c.s. of the reflection)
882
883 G4ThreeVector es1, es2, es3;
884
885 // z-Coordinate is the surface normal, has already length 1
886
887 es3 = Normal;
888
889 // Perpendicular part of incident direction w.r.t. normal
890 es1 = direction.perpPart(Normal);
891
892 // Set to unit length: x-Coordinate
893
894 es1.setMag(1.);
895 es2 = es1;
896
897 // y-Coordinate is the pi/2-rotation of x-axis around z-axis
898
899 es2.rotate(90.*degree, es3);
900
901 // Transformation matrix consists just of the three coordinate vectors
902 // as the global coordinate system is the 'standard' coordinate system
903
904 G4RotationMatrix matrix(es1, es2, es3);
905
906 return matrix;
907 }
908
909 void G4UCNBoundaryProcess::BoundaryProcessVerbose() const
910 {
911 if ( theStatus == Undefined )
912 G4cout << " *** Undefined *** " << G4endl;
913 if ( theStatus == NotAtBoundary )
914 G4cout << " *** NotAtBoundary *** " << G4endl;
915 if ( theStatus == SameMaterial )
916 G4cout << " *** SameMaterial *** " << G4endl;
917 if ( theStatus == StepTooSmall )
918 G4cout << " *** StepTooSmall *** " << G4endl;
919 if ( theStatus == NoMPT )
920 G4cout << " *** No G4UCNMaterialPropertiesTable *** " << G4endl;
921 if ( theStatus == NoMRT )
922 G4cout << " *** No MicroRoughness Table *** " << G4endl;
923 if ( theStatus == NoMRCondition )
924 G4cout << " *** MicroRoughness Condition not satisfied *** " << G4endl;
925 if ( theStatus == Absorption )
926 G4cout << " *** Loss on Surface *** " << G4endl;
927 if ( theStatus == Ezero )
928 G4cout << " *** Ezero on Surface *** " << G4endl;
929 if ( theStatus == Flip )
930 G4cout << " *** Spin Flip on Surface *** " << G4endl;
931 if ( theStatus == SpecularReflection )
932 G4cout << " *** Specular Reflection *** " << G4endl;
933 if ( theStatus == LambertianReflection )
934 G4cout << " *** LambertianR Reflection *** " << G4endl;
935 if ( theStatus == MRDiffuseReflection )
936 G4cout << " *** MR Model Diffuse Reflection *** " << G4endl;
937 if ( theStatus == SnellTransmit )
938 G4cout << " *** Snell Transmission *** " << G4endl;
939 if ( theStatus == MRDiffuseTransmit )
940 G4cout << " *** MR Model Diffuse Transmission *** " << G4endl;
941 }
942
943 void G4UCNBoundaryProcess::BoundaryProcessSummary(void) const
944 {
945 G4cout << "Sum NoMT: "
946 << nNoMPT << G4endl;
947 G4cout << "Sum NoMRT: "
948 << nNoMRT << G4endl;
949 G4cout << "Sum NoMRCondition: "
950 << nNoMRCondition << G4endl;
951 G4cout << "Sum No. E < V Loss: "
952 << nAbsorption << G4endl;
953 G4cout << "Sum No. E > V Ezero: "
954 << nEzero << G4endl;
955 G4cout << "Sum No. E < V SpinFlip: "
956 << nFlip << G4endl;
957 G4cout << "Sum No. E > V Specular Reflection: "
958 << aSpecularReflection << G4endl;
959 G4cout << "Sum No. E < V Specular Reflection: "
960 << bSpecularReflection << G4endl;
961 G4cout << "Sum No. E < V Lambertian Reflection: "
962 << bLambertianReflection << G4endl;
963 G4cout << "Sum No. E > V MR DiffuseReflection: "
964 << aMRDiffuseReflection << G4endl;
965 G4cout << "Sum No. E < V MR DiffuseReflection: "
966 << bMRDiffuseReflection << G4endl;
967 G4cout << "Sum No. E > V SnellTransmit: "
968 << nSnellTransmit << G4endl;
969 G4cout << "Sum No. E > V MR SnellTransmit: "
970 << mSnellTransmit << G4endl;
971 G4cout << "Sum No. E > V DiffuseTransmit: "
972 << aMRDiffuseTransmit << G4endl;
973 G4cout << " " << G4endl;
974 }
975
976 G4bool G4UCNBoundaryProcess::InvokeSD(const G4Step* pStep)
977 {
978 G4Step aStep = *pStep;
979
980 aStep.AddTotalEnergyDeposit(pStep->GetTrack()->GetKineticEnergy());
981
982 G4VSensitiveDetector* sd = aStep.GetPostStepPoint()->GetSensitiveDetector();
983 if (sd) return sd->Hit(&aStep);
984 else return false;
985 }
986