Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4OpBoundaryProcess.cc
이 파일의 문서화 페이지로 가기
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 //
27 // Optical Photon Boundary Process Class Implementation
29 //
30 // File: G4OpBoundaryProcess.cc
31 // Description: Discrete Process -- reflection/refraction at
32 // optical interfaces
33 // Version: 1.1
34 // Created: 1997-06-18
35 // Modified: 1998-05-25 - Correct parallel component of polarization
36 // (thanks to: Stefano Magni + Giovanni Pieri)
37 // 1998-05-28 - NULL Rindex pointer before reuse
38 // (thanks to: Stefano Magni)
39 // 1998-06-11 - delete *sint1 in oblique reflection
40 // (thanks to: Giovanni Pieri)
41 // 1998-06-19 - move from GetLocalExitNormal() to the new
42 // method: GetLocalExitNormal(&valid) to get
43 // the surface normal in all cases
44 // 1998-11-07 - NULL OpticalSurface pointer before use
45 // comparison not sharp for: std::abs(cost1) < 1.0
46 // remove sin1, sin2 in lines 556,567
47 // (thanks to Stefano Magni)
48 // 1999-10-10 - Accommodate changes done in DoAbsorption by
49 // changing logic in DielectricMetal
50 // 2001-10-18 - avoid Linux (gcc-2.95.2) warning about variables
51 // might be used uninitialized in this function
52 // moved E2_perp, E2_parl and E2_total out of 'if'
53 // 2003-11-27 - Modified line 168-9 to reflect changes made to
54 // G4OpticalSurface class ( by Fan Lei)
55 // 2004-02-02 - Set theStatus = Undefined at start of DoIt
56 // 2005-07-28 - add G4ProcessType to constructor
57 // 2006-11-04 - add capability of calculating the reflectivity
58 // off a metal surface by way of a complex index
59 // of refraction - Thanks to Sehwook Lee and John
60 // Hauptman (Dept. of Physics - Iowa State Univ.)
61 // 2009-11-10 - add capability of simulating surface reflections
62 // with Look-Up-Tables (LUT) containing measured
63 // optical reflectance for a variety of surface
64 // treatments - Thanks to Martin Janecek and
65 // William Moses (Lawrence Berkeley National Lab.)
66 // 2013-06-01 - add the capability of simulating the transmission
67 // of a dichronic filter
68 // 2017-02-24 - add capability of simulating surface reflections
69 // with Look-Up-Tables (LUT) developed in DAVIS
70 //
71 // Author: Peter Gumplinger
72 // adopted from work by Werner Keil - April 2/96
73 // mail: gum@triumf.ca
74 //
76 
77 #include "G4ios.hh"
78 #include "G4PhysicalConstants.hh"
79 #include "G4OpProcessSubType.hh"
80 
81 #include "G4OpBoundaryProcess.hh"
82 #include "G4GeometryTolerance.hh"
83 
84 #include "G4VSensitiveDetector.hh"
85 #include "G4ParallelWorldProcess.hh"
86 
87 #include "G4SystemOfUnits.hh"
88 
90 // Class Implementation
92 
94  // Operators
96 
97 // G4OpBoundaryProcess::operator=(const G4OpBoundaryProcess &right)
98 // {
99 // }
100 
102  // Constructors
104 
105 G4OpBoundaryProcess::G4OpBoundaryProcess(const G4String& processName,
106  G4ProcessType type)
107  : G4VDiscreteProcess(processName, type)
108 {
109  if ( verboseLevel > 0) {
110  G4cout << GetProcessName() << " is created " << G4endl;
111  }
112 
113  SetProcessSubType(fOpBoundary);
114 
115  theStatus = Undefined;
116  theModel = glisur;
117  theFinish = polished;
118  theReflectivity = 1.;
119  theEfficiency = 0.;
120  theTransmittance = 0.;
121 
122  theSurfaceRoughness = 0.;
123 
124  prob_sl = 0.;
125  prob_ss = 0.;
126  prob_bs = 0.;
127 
128  PropertyPointer = NULL;
129  PropertyPointer1 = NULL;
130  PropertyPointer2 = NULL;
131 
132  Material1 = NULL;
133  Material2 = NULL;
134 
135  OpticalSurface = NULL;
136 
137  kCarTolerance = G4GeometryTolerance::GetInstance()
138  ->GetSurfaceTolerance();
139 
140  iTE = iTM = 0;
141  thePhotonMomentum = 0.;
142  Rindex1 = Rindex2 = 1.;
143  cost1 = cost2 = sint1 = sint2 = 0.;
144 
145  idx = idy = 0;
146  DichroicVector = NULL;
147 
148  fInvokeSD = true;
149 }
150 
151 // G4OpBoundaryProcess::G4OpBoundaryProcess(const G4OpBoundaryProcess &right)
152 // {
153 // }
154 
156  // Destructors
158 
159 G4OpBoundaryProcess::~G4OpBoundaryProcess(){}
160 
162  // Methods
164 
165 // PostStepDoIt
166 // ------------
167 //
168 
169 G4VParticleChange*
170 G4OpBoundaryProcess::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
171 {
172  theStatus = Undefined;
173 
174  aParticleChange.Initialize(aTrack);
175  aParticleChange.ProposeVelocity(aTrack.GetVelocity());
176 
177  // Get hyperStep from G4ParallelWorldProcess
178  // NOTE: PostSetpDoIt of this process should be
179  // invoked after G4ParallelWorldProcess!
180 
181  const G4Step* pStep = &aStep;
182 
183  const G4Step* hStep = G4ParallelWorldProcess::GetHyperStep();
184 
185  if (hStep) pStep = hStep;
186 
187  G4bool isOnBoundary =
188  (pStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary);
189 
190  if (isOnBoundary) {
191  Material1 = pStep->GetPreStepPoint()->GetMaterial();
192  Material2 = pStep->GetPostStepPoint()->GetMaterial();
193  } else {
194  theStatus = NotAtBoundary;
195  if ( verboseLevel > 0) BoundaryProcessVerbose();
196  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
197  }
198 
199  G4VPhysicalVolume* thePrePV =
200  pStep->GetPreStepPoint() ->GetPhysicalVolume();
201  G4VPhysicalVolume* thePostPV =
202  pStep->GetPostStepPoint()->GetPhysicalVolume();
203 
204  if ( verboseLevel > 0 ) {
205  G4cout << " Photon at Boundary! " << G4endl;
206  if (thePrePV) G4cout << " thePrePV: " << thePrePV->GetName() << G4endl;
207  if (thePostPV) G4cout << " thePostPV: " << thePostPV->GetName() << G4endl;
208  }
209 
210  if (aTrack.GetStepLength()<=kCarTolerance/2){
211  theStatus = StepTooSmall;
212  if ( verboseLevel > 0) BoundaryProcessVerbose();
213  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
214  }
215 
216  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
217 
218  thePhotonMomentum = aParticle->GetTotalMomentum();
219  OldMomentum = aParticle->GetMomentumDirection();
220  OldPolarization = aParticle->GetPolarization();
221 
222  if ( verboseLevel > 0 ) {
223  G4cout << " Old Momentum Direction: " << OldMomentum << G4endl;
224  G4cout << " Old Polarization: " << OldPolarization << G4endl;
225  }
226 
227  G4ThreeVector theGlobalPoint = pStep->GetPostStepPoint()->GetPosition();
228 
229  G4bool valid;
230  // Use the new method for Exit Normal in global coordinates,
231  // which provides the normal more reliably.
232 
233  // ID of Navigator which limits step
234 
235  G4int hNavId = G4ParallelWorldProcess::GetHypNavigatorID();
236  std::vector<G4Navigator*>::iterator iNav =
237  G4TransportationManager::GetTransportationManager()->
238  GetActiveNavigatorsIterator();
239  theGlobalNormal =
240  (iNav[hNavId])->GetGlobalExitNormal(theGlobalPoint,&valid);
241 
242  if (valid) {
243  theGlobalNormal = -theGlobalNormal;
244  }
245  else
246  {
247  G4ExceptionDescription ed;
248  ed << " G4OpBoundaryProcess/PostStepDoIt(): "
249  << " The Navigator reports that it returned an invalid normal"
250  << G4endl;
251  G4Exception("G4OpBoundaryProcess::PostStepDoIt", "OpBoun01",
252  EventMustBeAborted,ed,
253  "Invalid Surface Normal - Geometry must return valid surface normal");
254  }
255 
256  if (OldMomentum * theGlobalNormal > 0.0) {
257 #ifdef G4OPTICAL_DEBUG
258  G4ExceptionDescription ed;
259  ed << " G4OpBoundaryProcess/PostStepDoIt(): "
260  << " theGlobalNormal points in a wrong direction. "
261  << G4endl;
262  ed << " The momentum of the photon arriving at interface (oldMomentum)"
263  << " must exit the volume cross in the step. " << G4endl;
264  ed << " So it MUST have dot < 0 with the normal that Exits the new volume (globalNormal)." << G4endl;
265  ed << " >> The dot product of oldMomentum and global Normal is " << OldMomentum*theGlobalNormal << G4endl;
266  ed << " Old Momentum (during step) = " << OldMomentum << G4endl;
267  ed << " Global Normal (Exiting New Vol) = " << theGlobalNormal << G4endl;
268  ed << G4endl;
269  G4Exception("G4OpBoundaryProcess::PostStepDoIt", "OpBoun02",
270  EventMustBeAborted, // Or JustWarning to see if it happens repeatedbly on one ray
271  ed,
272  "Invalid Surface Normal - Geometry must return valid surface normal pointing in the right direction");
273 #else
274  theGlobalNormal = -theGlobalNormal;
275 #endif
276  }
277 
278  G4MaterialPropertiesTable* aMaterialPropertiesTable;
279  G4MaterialPropertyVector* Rindex;
280 
281  aMaterialPropertiesTable = Material1->GetMaterialPropertiesTable();
282  if (aMaterialPropertiesTable) {
283  Rindex = aMaterialPropertiesTable->GetProperty(kRINDEX);
284  }
285  else {
286  theStatus = NoRINDEX;
287  if ( verboseLevel > 0) BoundaryProcessVerbose();
288  aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
289  aParticleChange.ProposeTrackStatus(fStopAndKill);
290  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
291  }
292 
293  if (Rindex) {
294  Rindex1 = Rindex->Value(thePhotonMomentum);
295  }
296  else {
297  theStatus = NoRINDEX;
298  if ( verboseLevel > 0) BoundaryProcessVerbose();
299  aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
300  aParticleChange.ProposeTrackStatus(fStopAndKill);
301  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
302  }
303 
304  theReflectivity = 1.;
305  theEfficiency = 0.;
306  theTransmittance = 0.;
307 
308  theSurfaceRoughness = 0.;
309 
310  theModel = glisur;
311  theFinish = polished;
312 
313  G4SurfaceType type = dielectric_dielectric;
314 
315  Rindex = NULL;
316  OpticalSurface = NULL;
317 
318  G4LogicalSurface* Surface = NULL;
319 
320  Surface = G4LogicalBorderSurface::GetSurface(thePrePV, thePostPV);
321 
322  if (Surface == NULL){
323  G4bool enteredDaughter= (thePostPV->GetMotherLogical() ==
324  thePrePV ->GetLogicalVolume());
325  if(enteredDaughter){
326  Surface =
327  G4LogicalSkinSurface::GetSurface(thePostPV->GetLogicalVolume());
328  if(Surface == NULL)
329  Surface =
330  G4LogicalSkinSurface::GetSurface(thePrePV->GetLogicalVolume());
331  }
332  else {
333  Surface =
334  G4LogicalSkinSurface::GetSurface(thePrePV->GetLogicalVolume());
335  if(Surface == NULL)
336  Surface =
337  G4LogicalSkinSurface::GetSurface(thePostPV->GetLogicalVolume());
338  }
339  }
340 
341  if (Surface) OpticalSurface =
342  dynamic_cast <G4OpticalSurface*> (Surface->GetSurfaceProperty());
343 
344  if (OpticalSurface) {
345 
346  type = OpticalSurface->GetType();
347  theModel = OpticalSurface->GetModel();
348  theFinish = OpticalSurface->GetFinish();
349 
350  aMaterialPropertiesTable = OpticalSurface->
351  GetMaterialPropertiesTable();
352 
353  if (aMaterialPropertiesTable) {
354 
355  if (theFinish == polishedbackpainted ||
356  theFinish == groundbackpainted ) {
357  Rindex = aMaterialPropertiesTable->GetProperty(kRINDEX);
358  if (Rindex) {
359  Rindex2 = Rindex->Value(thePhotonMomentum);
360  }
361  else {
362  theStatus = NoRINDEX;
363  if ( verboseLevel > 0) BoundaryProcessVerbose();
364  aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
365  aParticleChange.ProposeTrackStatus(fStopAndKill);
366  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
367  }
368  }
369 
370  PropertyPointer =
371  aMaterialPropertiesTable->GetProperty(kREFLECTIVITY);
372  PropertyPointer1 =
373  aMaterialPropertiesTable->GetProperty(kREALRINDEX);
374  PropertyPointer2 =
375  aMaterialPropertiesTable->GetProperty(kIMAGINARYRINDEX);
376 
377  iTE = 1;
378  iTM = 1;
379 
380  if (PropertyPointer) {
381 
382  theReflectivity =
383  PropertyPointer->Value(thePhotonMomentum);
384 
385  } else if (PropertyPointer1 && PropertyPointer2) {
386 
387  CalculateReflectivity();
388 
389  }
390 
391  PropertyPointer =
392  aMaterialPropertiesTable->GetProperty(kEFFICIENCY);
393  if (PropertyPointer) {
394  theEfficiency =
395  PropertyPointer->Value(thePhotonMomentum);
396  }
397 
398  PropertyPointer =
399  aMaterialPropertiesTable->GetProperty(kTRANSMITTANCE);
400  if (PropertyPointer) {
401  theTransmittance =
402  PropertyPointer->Value(thePhotonMomentum);
403  }
404 
405  if (aMaterialPropertiesTable->
406  ConstPropertyExists("SURFACEROUGHNESS"))
407  theSurfaceRoughness = aMaterialPropertiesTable->
408  GetConstProperty(kSURFACEROUGHNESS);
409 
410  if ( theModel == unified ) {
411  PropertyPointer =
412  aMaterialPropertiesTable->GetProperty(kSPECULARLOBECONSTANT);
413  if (PropertyPointer) {
414  prob_sl =
415  PropertyPointer->Value(thePhotonMomentum);
416  } else {
417  prob_sl = 0.0;
418  }
419 
420  PropertyPointer =
421  aMaterialPropertiesTable->GetProperty(kSPECULARSPIKECONSTANT);
422  if (PropertyPointer) {
423  prob_ss =
424  PropertyPointer->Value(thePhotonMomentum);
425  } else {
426  prob_ss = 0.0;
427  }
428 
429  PropertyPointer =
430  aMaterialPropertiesTable->GetProperty(kBACKSCATTERCONSTANT);
431  if (PropertyPointer) {
432  prob_bs =
433  PropertyPointer->Value(thePhotonMomentum);
434  } else {
435  prob_bs = 0.0;
436  }
437  }
438  }
439  else if (theFinish == polishedbackpainted ||
440  theFinish == groundbackpainted ) {
441  aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
442  aParticleChange.ProposeTrackStatus(fStopAndKill);
443  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
444  }
445  }
446 
447  if (type == dielectric_dielectric ) {
448  if (theFinish == polished || theFinish == ground ) {
449 
450  if (Material1 == Material2){
451  theStatus = SameMaterial;
452  if ( verboseLevel > 0) BoundaryProcessVerbose();
453  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
454  }
455  aMaterialPropertiesTable =
456  Material2->GetMaterialPropertiesTable();
457  if (aMaterialPropertiesTable)
458  Rindex = aMaterialPropertiesTable->GetProperty(kRINDEX);
459  if (Rindex) {
460  Rindex2 = Rindex->Value(thePhotonMomentum);
461  }
462  else {
463  theStatus = NoRINDEX;
464  if ( verboseLevel > 0) BoundaryProcessVerbose();
465  aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
466  aParticleChange.ProposeTrackStatus(fStopAndKill);
467  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
468  }
469  }
470  }
471 
472  if (type == dielectric_metal) {
473 
474  DielectricMetal();
475 
476  }
477  else if (type == dielectric_LUT) {
478 
479  DielectricLUT();
480 
481  }
482  else if (type == dielectric_LUTDAVIS) {
483 
484  DielectricLUTDAVIS();
485 
486  }
487  else if (type == dielectric_dichroic) {
488 
489  DielectricDichroic();
490 
491  }
492  else if (type == dielectric_dielectric) {
493 
494  if ( theFinish == polishedbackpainted ||
495  theFinish == groundbackpainted ) {
496  DielectricDielectric();
497  }
498  else {
499  G4double rand = G4UniformRand();
500  if ( rand > theReflectivity ) {
501  if (rand > theReflectivity + theTransmittance) {
502  DoAbsorption();
503  } else {
504  theStatus = Transmission;
505  NewMomentum = OldMomentum;
506  NewPolarization = OldPolarization;
507  }
508  }
509  else {
510  if ( theFinish == polishedfrontpainted ) {
511  DoReflection();
512  }
513  else if ( theFinish == groundfrontpainted ) {
514  theStatus = LambertianReflection;
515  DoReflection();
516  }
517  else {
518  DielectricDielectric();
519  }
520  }
521  }
522  }
523  else {
524 
525  G4cerr << " Error: G4BoundaryProcess: illegal boundary type " << G4endl;
526  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
527 
528  }
529 
530  NewMomentum = NewMomentum.unit();
531  NewPolarization = NewPolarization.unit();
532 
533  if ( verboseLevel > 0) {
534  G4cout << " New Momentum Direction: " << NewMomentum << G4endl;
535  G4cout << " New Polarization: " << NewPolarization << G4endl;
536  BoundaryProcessVerbose();
537  }
538 
539  aParticleChange.ProposeMomentumDirection(NewMomentum);
540  aParticleChange.ProposePolarization(NewPolarization);
541 
542  if ( theStatus == FresnelRefraction || theStatus == Transmission ) {
543  G4MaterialPropertyVector* groupvel =
544  Material2->GetMaterialPropertiesTable()->GetProperty(kGROUPVEL);
545  G4double finalVelocity = groupvel->Value(thePhotonMomentum);
546  aParticleChange.ProposeVelocity(finalVelocity);
547  }
548 
549  if ( theStatus == Detection && fInvokeSD ) InvokeSD(pStep);
550 
551  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
552 }
553 
554 void G4OpBoundaryProcess::BoundaryProcessVerbose() const
555 {
556  if ( theStatus == Undefined )
557  G4cout << " *** Undefined *** " << G4endl;
558  if ( theStatus == Transmission )
559  G4cout << " *** Transmission *** " << G4endl;
560  if ( theStatus == FresnelRefraction )
561  G4cout << " *** FresnelRefraction *** " << G4endl;
562  if ( theStatus == FresnelReflection )
563  G4cout << " *** FresnelReflection *** " << G4endl;
564  if ( theStatus == TotalInternalReflection )
565  G4cout << " *** TotalInternalReflection *** " << G4endl;
566  if ( theStatus == LambertianReflection )
567  G4cout << " *** LambertianReflection *** " << G4endl;
568  if ( theStatus == LobeReflection )
569  G4cout << " *** LobeReflection *** " << G4endl;
570  if ( theStatus == SpikeReflection )
571  G4cout << " *** SpikeReflection *** " << G4endl;
572  if ( theStatus == BackScattering )
573  G4cout << " *** BackScattering *** " << G4endl;
574  if ( theStatus == PolishedLumirrorAirReflection )
575  G4cout << " *** PolishedLumirrorAirReflection *** " << G4endl;
576  if ( theStatus == PolishedLumirrorGlueReflection )
577  G4cout << " *** PolishedLumirrorGlueReflection *** " << G4endl;
578  if ( theStatus == PolishedAirReflection )
579  G4cout << " *** PolishedAirReflection *** " << G4endl;
580  if ( theStatus == PolishedTeflonAirReflection )
581  G4cout << " *** PolishedTeflonAirReflection *** " << G4endl;
582  if ( theStatus == PolishedTiOAirReflection )
583  G4cout << " *** PolishedTiOAirReflection *** " << G4endl;
584  if ( theStatus == PolishedTyvekAirReflection )
585  G4cout << " *** PolishedTyvekAirReflection *** " << G4endl;
586  if ( theStatus == PolishedVM2000AirReflection )
587  G4cout << " *** PolishedVM2000AirReflection *** " << G4endl;
588  if ( theStatus == PolishedVM2000GlueReflection )
589  G4cout << " *** PolishedVM2000GlueReflection *** " << G4endl;
590  if ( theStatus == EtchedLumirrorAirReflection )
591  G4cout << " *** EtchedLumirrorAirReflection *** " << G4endl;
592  if ( theStatus == EtchedLumirrorGlueReflection )
593  G4cout << " *** EtchedLumirrorGlueReflection *** " << G4endl;
594  if ( theStatus == EtchedAirReflection )
595  G4cout << " *** EtchedAirReflection *** " << G4endl;
596  if ( theStatus == EtchedTeflonAirReflection )
597  G4cout << " *** EtchedTeflonAirReflection *** " << G4endl;
598  if ( theStatus == EtchedTiOAirReflection )
599  G4cout << " *** EtchedTiOAirReflection *** " << G4endl;
600  if ( theStatus == EtchedTyvekAirReflection )
601  G4cout << " *** EtchedTyvekAirReflection *** " << G4endl;
602  if ( theStatus == EtchedVM2000AirReflection )
603  G4cout << " *** EtchedVM2000AirReflection *** " << G4endl;
604  if ( theStatus == EtchedVM2000GlueReflection )
605  G4cout << " *** EtchedVM2000GlueReflection *** " << G4endl;
606  if ( theStatus == GroundLumirrorAirReflection )
607  G4cout << " *** GroundLumirrorAirReflection *** " << G4endl;
608  if ( theStatus == GroundLumirrorGlueReflection )
609  G4cout << " *** GroundLumirrorGlueReflection *** " << G4endl;
610  if ( theStatus == GroundAirReflection )
611  G4cout << " *** GroundAirReflection *** " << G4endl;
612  if ( theStatus == GroundTeflonAirReflection )
613  G4cout << " *** GroundTeflonAirReflection *** " << G4endl;
614  if ( theStatus == GroundTiOAirReflection )
615  G4cout << " *** GroundTiOAirReflection *** " << G4endl;
616  if ( theStatus == GroundTyvekAirReflection )
617  G4cout << " *** GroundTyvekAirReflection *** " << G4endl;
618  if ( theStatus == GroundVM2000AirReflection )
619  G4cout << " *** GroundVM2000AirReflection *** " << G4endl;
620  if ( theStatus == GroundVM2000GlueReflection )
621  G4cout << " *** GroundVM2000GlueReflection *** " << G4endl;
622  if ( theStatus == Absorption )
623  G4cout << " *** Absorption *** " << G4endl;
624  if ( theStatus == Detection )
625  G4cout << " *** Detection *** " << G4endl;
626  if ( theStatus == NotAtBoundary )
627  G4cout << " *** NotAtBoundary *** " << G4endl;
628  if ( theStatus == SameMaterial )
629  G4cout << " *** SameMaterial *** " << G4endl;
630  if ( theStatus == StepTooSmall )
631  G4cout << " *** StepTooSmall *** " << G4endl;
632  if ( theStatus == NoRINDEX )
633  G4cout << " *** NoRINDEX *** " << G4endl;
634  if ( theStatus == Dichroic )
635  G4cout << " *** Dichroic Transmission *** " << G4endl;
636 }
637 
638 G4ThreeVector
639 G4OpBoundaryProcess::GetFacetNormal(const G4ThreeVector& Momentum,
640  const G4ThreeVector& Normal ) const
641 {
642  G4ThreeVector FacetNormal;
643 
644  if (theModel == unified || theModel == LUT || theModel== DAVIS) {
645 
646  /* This function code alpha to a random value taken from the
647  distribution p(alpha) = g(alpha; 0, sigma_alpha)*std::sin(alpha),
648  for alpha > 0 and alpha < 90, where g(alpha; 0, sigma_alpha)
649  is a gaussian distribution with mean 0 and standard deviation
650  sigma_alpha. */
651 
652  G4double alpha;
653 
654  G4double sigma_alpha = 0.0;
655  if (OpticalSurface) sigma_alpha = OpticalSurface->GetSigmaAlpha();
656 
657  if (sigma_alpha == 0.0) return FacetNormal = Normal;
658 
659  G4double f_max = std::min(1.0,4.*sigma_alpha);
660 
661  G4double phi, SinAlpha, CosAlpha, SinPhi, CosPhi, unit_x, unit_y, unit_z;
662  G4ThreeVector tmpNormal;
663 
664  do {
665  do {
666  alpha = G4RandGauss::shoot(0.0,sigma_alpha);
667  // Loop checking, 13-Aug-2015, Peter Gumplinger
668  } while (G4UniformRand()*f_max > std::sin(alpha) || alpha >= halfpi );
669 
670  phi = G4UniformRand()*twopi;
671 
672  SinAlpha = std::sin(alpha);
673  CosAlpha = std::cos(alpha);
674  SinPhi = std::sin(phi);
675  CosPhi = std::cos(phi);
676 
677  unit_x = SinAlpha * CosPhi;
678  unit_y = SinAlpha * SinPhi;
679  unit_z = CosAlpha;
680 
681  FacetNormal.setX(unit_x);
682  FacetNormal.setY(unit_y);
683  FacetNormal.setZ(unit_z);
684 
685  tmpNormal = Normal;
686 
687  FacetNormal.rotateUz(tmpNormal);
688  // Loop checking, 13-Aug-2015, Peter Gumplinger
689  } while (Momentum * FacetNormal >= 0.0);
690  }
691  else {
692 
693  G4double polish = 1.0;
694  if (OpticalSurface) polish = OpticalSurface->GetPolish();
695 
696  if (polish < 1.0) {
697  do {
698  G4ThreeVector smear;
699  do {
700  smear.setX(2.*G4UniformRand()-1.0);
701  smear.setY(2.*G4UniformRand()-1.0);
702  smear.setZ(2.*G4UniformRand()-1.0);
703  // Loop checking, 13-Aug-2015, Peter Gumplinger
704  } while (smear.mag()>1.0);
705  smear = (1.-polish) * smear;
706  FacetNormal = Normal + smear;
707  // Loop checking, 13-Aug-2015, Peter Gumplinger
708  } while (Momentum * FacetNormal >= 0.0);
709  FacetNormal = FacetNormal.unit();
710  }
711  else {
712  FacetNormal = Normal;
713  }
714  }
715  return FacetNormal;
716 }
717 
718 void G4OpBoundaryProcess::DielectricMetal()
719 {
720  G4int n = 0;
721  G4double rand, PdotN, EdotN;
722  G4ThreeVector A_trans, A_paral;
723 
724  do {
725 
726  n++;
727 
728  rand = G4UniformRand();
729  if ( rand > theReflectivity && n == 1 ) {
730  if (rand > theReflectivity + theTransmittance) {
731  DoAbsorption();
732  } else {
733  theStatus = Transmission;
734  NewMomentum = OldMomentum;
735  NewPolarization = OldPolarization;
736  }
737  break;
738  }
739  else {
740 
741  if (PropertyPointer1 && PropertyPointer2) {
742  if ( n > 1 ) {
743  CalculateReflectivity();
744  if ( !G4BooleanRand(theReflectivity) ) {
745  DoAbsorption();
746  break;
747  }
748  }
749  }
750 
751  if ( theModel == glisur || theFinish == polished ) {
752 
753  DoReflection();
754 
755  } else {
756 
757  if ( n == 1 ) ChooseReflection();
758 
759  if ( theStatus == LambertianReflection ) {
760  DoReflection();
761  }
762  else if ( theStatus == BackScattering ) {
763  NewMomentum = -OldMomentum;
764  NewPolarization = -OldPolarization;
765  }
766  else {
767 
768  if(theStatus==LobeReflection){
769  if ( PropertyPointer1 && PropertyPointer2 ){
770  } else {
771  theFacetNormal =
772  GetFacetNormal(OldMomentum,theGlobalNormal);
773  }
774  }
775 
776  PdotN = OldMomentum * theFacetNormal;
777  NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
778  EdotN = OldPolarization * theFacetNormal;
779 
780  if (sint1 > 0.0 ) {
781  A_trans = OldMomentum.cross(theFacetNormal);
782  A_trans = A_trans.unit();
783  } else {
784  A_trans = OldPolarization;
785  }
786  A_paral = NewMomentum.cross(A_trans);
787  A_paral = A_paral.unit();
788 
789  if(iTE>0&&iTM>0) {
790  NewPolarization =
791  -OldPolarization + (2.*EdotN)*theFacetNormal;
792  } else if (iTE>0) {
793  NewPolarization = -A_trans;
794  } else if (iTM>0) {
795  NewPolarization = -A_paral;
796  }
797 
798  }
799 
800  }
801 
802  OldMomentum = NewMomentum;
803  OldPolarization = NewPolarization;
804 
805  }
806 
807  // Loop checking, 13-Aug-2015, Peter Gumplinger
808  } while (NewMomentum * theGlobalNormal < 0.0);
809 }
810 
811 void G4OpBoundaryProcess::DielectricLUT()
812 {
813  G4int thetaIndex, phiIndex;
814  G4double AngularDistributionValue, thetaRad, phiRad, EdotN;
815  G4ThreeVector PerpendicularVectorTheta, PerpendicularVectorPhi;
816 
817  theStatus = G4OpBoundaryProcessStatus(G4int(theFinish) +
818  (G4int(NoRINDEX)-G4int(groundbackpainted)));
819 
820  G4int thetaIndexMax = OpticalSurface->GetThetaIndexMax();
821  G4int phiIndexMax = OpticalSurface->GetPhiIndexMax();
822 
823  G4double rand;
824 
825  do {
826  rand = G4UniformRand();
827  if ( rand > theReflectivity ) {
828  if (rand > theReflectivity + theTransmittance) {
829  DoAbsorption();
830  } else {
831  theStatus = Transmission;
832  NewMomentum = OldMomentum;
833  NewPolarization = OldPolarization;
834  }
835  break;
836  }
837  else {
838  // Calculate Angle between Normal and Photon Momentum
839  G4double anglePhotonToNormal =
840  OldMomentum.angle(-theGlobalNormal);
841  // Round it to closest integer
842  G4int angleIncident = G4int(std::floor(180/pi*anglePhotonToNormal+0.5));
843 
844  // Take random angles THETA and PHI,
845  // and see if below Probability - if not - Redo
846  do {
847  thetaIndex = G4RandFlat::shootInt(thetaIndexMax-1);
848  phiIndex = G4RandFlat::shootInt(phiIndexMax-1);
849  // Find probability with the new indeces from LUT
850  AngularDistributionValue = OpticalSurface ->
851  GetAngularDistributionValue(angleIncident,
852  thetaIndex,
853  phiIndex);
854  // Loop checking, 13-Aug-2015, Peter Gumplinger
855  } while ( !G4BooleanRand(AngularDistributionValue) );
856 
857  thetaRad = (-90 + 4*thetaIndex)*pi/180;
858  phiRad = (-90 + 5*phiIndex)*pi/180;
859  // Rotate Photon Momentum in Theta, then in Phi
860  NewMomentum = -OldMomentum;
861 
862  PerpendicularVectorTheta = NewMomentum.cross(theGlobalNormal);
863  if (PerpendicularVectorTheta.mag() < kCarTolerance )
864  PerpendicularVectorTheta = NewMomentum.orthogonal();
865  NewMomentum =
866  NewMomentum.rotate(anglePhotonToNormal-thetaRad,
867  PerpendicularVectorTheta);
868  PerpendicularVectorPhi =
869  PerpendicularVectorTheta.cross(NewMomentum);
870  NewMomentum = NewMomentum.rotate(-phiRad,PerpendicularVectorPhi);
871 
872  // Rotate Polarization too:
873  theFacetNormal = (NewMomentum - OldMomentum).unit();
874  EdotN = OldPolarization * theFacetNormal;
875  NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
876  }
877  // Loop checking, 13-Aug-2015, Peter Gumplinger
878  } while (NewMomentum * theGlobalNormal <= 0.0);
879 }
880 
881 void G4OpBoundaryProcess::DielectricLUTDAVIS()
882 {
883  G4int angindex, random, angleIncident;
884  G4double ReflectivityValue, elevation, azimuth, EdotN;
885  G4double anglePhotonToNormal;
886 
887  G4int LUTbin = OpticalSurface->GetLUTbins();
888 
889  G4double rand = G4UniformRand();
890 
891  do {
892 
893  anglePhotonToNormal = OldMomentum.angle(-theGlobalNormal);
894  angleIncident = G4int(std::floor(180/pi*anglePhotonToNormal+0.5));
895 
896  ReflectivityValue = OpticalSurface -> GetReflectivityLUTValue(angleIncident);
897 
898  if ( rand > ReflectivityValue ) {
899 
900  if ( theEfficiency > 0 ) {
901  DoAbsorption();
902  break;
903  }
904  else {
905 
906  theStatus = Transmission;
907 
908  if (angleIncident <= 0.01) {
909  NewMomentum = OldMomentum;
910  break;
911 
912  }
913 
914  do {
915  random = G4RandFlat::shootInt(1,LUTbin+1);
916  angindex = (((random*2)-1))+angleIncident*LUTbin*2 + 3640000;
917 
918  azimuth = OpticalSurface -> GetAngularDistributionValueLUT(angindex-1);
919  elevation= OpticalSurface -> GetAngularDistributionValueLUT(angindex);
920 
921  } while ( elevation == 0 && azimuth == 0);
922 
923  NewMomentum = -OldMomentum;
924 
925  G4ThreeVector v = theGlobalNormal.cross(-NewMomentum);
926  G4ThreeVector vNorm = v/v.mag();
927  G4ThreeVector u = vNorm.cross(theGlobalNormal);
928 
929  u = u *= (std::sin(elevation) * std::cos(azimuth));
930  v = vNorm *= (std::sin(elevation) * std::sin(azimuth));
931  G4ThreeVector w = theGlobalNormal *= (std::cos(elevation));
932  NewMomentum = G4ThreeVector(u+v+w);
933 
934  // Rotate Polarization too:
935  theFacetNormal = (NewMomentum - OldMomentum).unit();
936  EdotN = OldPolarization * theFacetNormal;
937  NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
938  }
939  }
940  else {
941 
942  theStatus = LobeReflection;
943 
944  if (angleIncident == 0) {
945  NewMomentum = -OldMomentum;
946  break;
947  }
948 
949  do {
950  random = G4RandFlat::shootInt(1,LUTbin+1);
951  angindex = (((random*2)-1))+(angleIncident-1)*LUTbin*2;
952 
953  azimuth = OpticalSurface -> GetAngularDistributionValueLUT(angindex-1);
954  elevation = OpticalSurface -> GetAngularDistributionValueLUT(angindex);
955  } while (elevation == 0 && azimuth == 0);
956 
957  NewMomentum = -OldMomentum;
958 
959  G4ThreeVector v = theGlobalNormal.cross(-NewMomentum);
960  G4ThreeVector vNorm = v/v.mag();
961  G4ThreeVector u = vNorm.cross(theGlobalNormal);
962 
963  u = u *= (std::sin(elevation) * std::cos(azimuth));
964  v = vNorm *= (std::sin(elevation) * std::sin(azimuth));
965  G4ThreeVector w = theGlobalNormal*=(std::cos(elevation));
966 
967  NewMomentum = G4ThreeVector(u+v+w);
968 
969  // Rotate Polarization too: (needs revision)
970  NewPolarization = OldPolarization;
971  }
972  } while (NewMomentum * theGlobalNormal <= 0.0);
973 }
974 
975 void G4OpBoundaryProcess::DielectricDichroic()
976 {
977  // Calculate Angle between Normal and Photon Momentum
978  G4double anglePhotonToNormal = OldMomentum.angle(-theGlobalNormal);
979 
980  // Round it to closest integer
981  G4double angleIncident = std::floor(180/pi*anglePhotonToNormal+0.5);
982 
983  if (!DichroicVector) {
984  if (OpticalSurface) DichroicVector = OpticalSurface->GetDichroicVector();
985  }
986 
987 
988  if (DichroicVector) {
989  G4double wavelength = h_Planck*c_light/thePhotonMomentum;
990  theTransmittance =
991  DichroicVector->Value(wavelength/nm,angleIncident,idx,idy)*perCent;
992 // G4cout << "wavelength: " << std::floor(wavelength/nm)
993 // << "nm" << G4endl;
994 // G4cout << "Incident angle: " << angleIncident << "deg" << G4endl;
995 // G4cout << "Transmittance: "
996 // << std::floor(theTransmittance/perCent) << "%" << G4endl;
997  } else {
998  G4ExceptionDescription ed;
999  ed << " G4OpBoundaryProcess/DielectricDichroic(): "
1000  << " The dichroic surface has no G4Physics2DVector"
1001  << G4endl;
1002  G4Exception("G4OpBoundaryProcess::DielectricDichroic", "OpBoun03",
1003  FatalException,ed,
1004  "A dichroic surface must have an associated G4Physics2DVector");
1005  }
1006 
1007  if ( !G4BooleanRand(theTransmittance) ) { // Not transmitted, so reflect
1008 
1009  if ( theModel == glisur || theFinish == polished ) {
1010  DoReflection();
1011  } else {
1012  ChooseReflection();
1013  if ( theStatus == LambertianReflection ) {
1014  DoReflection();
1015  } else if ( theStatus == BackScattering ) {
1016  NewMomentum = -OldMomentum;
1017  NewPolarization = -OldPolarization;
1018  } else {
1019  G4double PdotN, EdotN;
1020  do {
1021  if (theStatus==LobeReflection)
1022  theFacetNormal = GetFacetNormal(OldMomentum,theGlobalNormal);
1023  PdotN = OldMomentum * theFacetNormal;
1024  NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
1025  // Loop checking, 13-Aug-2015, Peter Gumplinger
1026  } while (NewMomentum * theGlobalNormal <= 0.0);
1027  EdotN = OldPolarization * theFacetNormal;
1028  NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
1029  }
1030  }
1031 
1032  } else {
1033 
1034  theStatus = Dichroic;
1035  NewMomentum = OldMomentum;
1036  NewPolarization = OldPolarization;
1037 
1038  }
1039 }
1040 
1041 void G4OpBoundaryProcess::DielectricDielectric()
1042 {
1043  G4bool Inside = false;
1044  G4bool Swap = false;
1045 
1046  G4bool SurfaceRoughnessCriterionPass = 1;
1047  if (theSurfaceRoughness != 0. && Rindex1 > Rindex2) {
1048  G4double wavelength = h_Planck*c_light/thePhotonMomentum;
1049  G4double SurfaceRoughnessCriterion =
1050  std::exp(-std::pow((4*pi*theSurfaceRoughness*Rindex1*cost1/wavelength),2));
1051  SurfaceRoughnessCriterionPass =
1052  G4BooleanRand(SurfaceRoughnessCriterion);
1053  }
1054 
1055  leap:
1056 
1057  G4bool Through = false;
1058  G4bool Done = false;
1059 
1060  G4double PdotN, EdotN;
1061 
1062  G4ThreeVector A_trans, A_paral, E1pp, E1pl;
1063  G4double E1_perp, E1_parl;
1064  G4double s1, s2, E2_perp, E2_parl, E2_total, TransCoeff;
1065  G4double E2_abs, C_parl, C_perp;
1066  G4double alpha;
1067 
1068  do {
1069 
1070  if (Through) {
1071  Swap = !Swap;
1072  Through = false;
1073  theGlobalNormal = -theGlobalNormal;
1074  G4SwapPtr(Material1,Material2);
1075  G4SwapObj(&Rindex1,&Rindex2);
1076  }
1077 
1078  if ( theFinish == polished ) {
1079  theFacetNormal = theGlobalNormal;
1080  }
1081  else {
1082  theFacetNormal =
1083  GetFacetNormal(OldMomentum,theGlobalNormal);
1084  }
1085 
1086  PdotN = OldMomentum * theFacetNormal;
1087  EdotN = OldPolarization * theFacetNormal;
1088 
1089  cost1 = - PdotN;
1090  if (std::abs(cost1) < 1.0-kCarTolerance){
1091  sint1 = std::sqrt(1.-cost1*cost1);
1092  sint2 = sint1*Rindex1/Rindex2; // *** Snell's Law ***
1093  }
1094  else {
1095  sint1 = 0.0;
1096  sint2 = 0.0;
1097  }
1098 
1099  if (sint2 >= 1.0) {
1100 
1101  // Simulate total internal reflection
1102 
1103  if (Swap) Swap = !Swap;
1104 
1105  theStatus = TotalInternalReflection;
1106 
1107  if ( !SurfaceRoughnessCriterionPass ) theStatus =
1108  LambertianReflection;
1109 
1110  if ( theModel == unified && theFinish != polished )
1111  ChooseReflection();
1112 
1113  if ( theStatus == LambertianReflection ) {
1114  DoReflection();
1115  }
1116  else if ( theStatus == BackScattering ) {
1117  NewMomentum = -OldMomentum;
1118  NewPolarization = -OldPolarization;
1119  }
1120  else {
1121 
1122  PdotN = OldMomentum * theFacetNormal;
1123  NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
1124  EdotN = OldPolarization * theFacetNormal;
1125  NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
1126 
1127  }
1128  }
1129  else if (sint2 < 1.0) {
1130 
1131  // Calculate amplitude for transmission (Q = P x N)
1132 
1133  if (cost1 > 0.0) {
1134  cost2 = std::sqrt(1.-sint2*sint2);
1135  }
1136  else {
1137  cost2 = -std::sqrt(1.-sint2*sint2);
1138  }
1139 
1140  if (sint1 > 0.0) {
1141  A_trans = OldMomentum.cross(theFacetNormal);
1142  A_trans = A_trans.unit();
1143  E1_perp = OldPolarization * A_trans;
1144  E1pp = E1_perp * A_trans;
1145  E1pl = OldPolarization - E1pp;
1146  E1_parl = E1pl.mag();
1147  }
1148  else {
1149  A_trans = OldPolarization;
1150  // Here we Follow Jackson's conventions and we set the
1151  // parallel component = 1 in case of a ray perpendicular
1152  // to the surface
1153  E1_perp = 0.0;
1154  E1_parl = 1.0;
1155  }
1156 
1157  s1 = Rindex1*cost1;
1158  E2_perp = 2.*s1*E1_perp/(Rindex1*cost1+Rindex2*cost2);
1159  E2_parl = 2.*s1*E1_parl/(Rindex2*cost1+Rindex1*cost2);
1160  E2_total = E2_perp*E2_perp + E2_parl*E2_parl;
1161  s2 = Rindex2*cost2*E2_total;
1162 
1163  if (theTransmittance > 0) TransCoeff = theTransmittance;
1164  else if (cost1 != 0.0) TransCoeff = s2/s1;
1165  else TransCoeff = 0.0;
1166 
1167  if ( !G4BooleanRand(TransCoeff) ) {
1168 
1169  // Simulate reflection
1170 
1171  if (Swap) Swap = !Swap;
1172 
1173  theStatus = FresnelReflection;
1174 
1175  if ( !SurfaceRoughnessCriterionPass ) theStatus =
1176  LambertianReflection;
1177 
1178  if ( theModel == unified && theFinish != polished )
1179  ChooseReflection();
1180 
1181  if ( theStatus == LambertianReflection ) {
1182  DoReflection();
1183  }
1184  else if ( theStatus == BackScattering ) {
1185  NewMomentum = -OldMomentum;
1186  NewPolarization = -OldPolarization;
1187  }
1188  else {
1189 
1190  PdotN = OldMomentum * theFacetNormal;
1191  NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
1192 
1193  if (sint1 > 0.0) { // incident ray oblique
1194 
1195  E2_parl = Rindex2*E2_parl/Rindex1 - E1_parl;
1196  E2_perp = E2_perp - E1_perp;
1197  E2_total = E2_perp*E2_perp + E2_parl*E2_parl;
1198  A_paral = NewMomentum.cross(A_trans);
1199  A_paral = A_paral.unit();
1200  E2_abs = std::sqrt(E2_total);
1201  C_parl = E2_parl/E2_abs;
1202  C_perp = E2_perp/E2_abs;
1203 
1204  NewPolarization = C_parl*A_paral + C_perp*A_trans;
1205 
1206  }
1207 
1208  else { // incident ray perpendicular
1209 
1210  if (Rindex2 > Rindex1) {
1211  NewPolarization = - OldPolarization;
1212  }
1213  else {
1214  NewPolarization = OldPolarization;
1215  }
1216 
1217  }
1218  }
1219  }
1220  else { // photon gets transmitted
1221 
1222  // Simulate transmission/refraction
1223 
1224  Inside = !Inside;
1225  Through = true;
1226  theStatus = FresnelRefraction;
1227 
1228  if (sint1 > 0.0) { // incident ray oblique
1229 
1230  alpha = cost1 - cost2*(Rindex2/Rindex1);
1231  NewMomentum = OldMomentum + alpha*theFacetNormal;
1232  NewMomentum = NewMomentum.unit();
1233 // PdotN = -cost2;
1234  A_paral = NewMomentum.cross(A_trans);
1235  A_paral = A_paral.unit();
1236  E2_abs = std::sqrt(E2_total);
1237  C_parl = E2_parl/E2_abs;
1238  C_perp = E2_perp/E2_abs;
1239 
1240  NewPolarization = C_parl*A_paral + C_perp*A_trans;
1241 
1242  }
1243  else { // incident ray perpendicular
1244 
1245  NewMomentum = OldMomentum;
1246  NewPolarization = OldPolarization;
1247 
1248  }
1249  }
1250  }
1251 
1252  OldMomentum = NewMomentum.unit();
1253  OldPolarization = NewPolarization.unit();
1254 
1255  if (theStatus == FresnelRefraction) {
1256  Done = (NewMomentum * theGlobalNormal <= 0.0);
1257  }
1258  else {
1259  Done = (NewMomentum * theGlobalNormal >= -kCarTolerance);
1260  }
1261 
1262  // Loop checking, 13-Aug-2015, Peter Gumplinger
1263  } while (!Done);
1264 
1265  if (Inside && !Swap) {
1266  if( theFinish == polishedbackpainted ||
1267  theFinish == groundbackpainted ) {
1268 
1269  G4double rand = G4UniformRand();
1270  if ( rand > theReflectivity ) {
1271  if (rand > theReflectivity + theTransmittance) {
1272  DoAbsorption();
1273  } else {
1274  theStatus = Transmission;
1275  NewMomentum = OldMomentum;
1276  NewPolarization = OldPolarization;
1277  }
1278  }
1279  else {
1280  if (theStatus != FresnelRefraction ) {
1281  theGlobalNormal = -theGlobalNormal;
1282  }
1283  else {
1284  Swap = !Swap;
1285  G4SwapPtr(Material1,Material2);
1286  G4SwapObj(&Rindex1,&Rindex2);
1287  }
1288  if ( theFinish == groundbackpainted )
1289  theStatus = LambertianReflection;
1290 
1291  DoReflection();
1292 
1293  theGlobalNormal = -theGlobalNormal;
1294  OldMomentum = NewMomentum;
1295 
1296  goto leap;
1297  }
1298  }
1299  }
1300 }
1301 
1302 // GetMeanFreePath
1303 // ---------------
1304 //
1305 G4double G4OpBoundaryProcess::GetMeanFreePath(const G4Track& ,
1306  G4double ,
1307  G4ForceCondition* condition)
1308 {
1309  *condition = Forced;
1310 
1311  return DBL_MAX;
1312 }
1313 
1314 G4double G4OpBoundaryProcess::GetIncidentAngle()
1315 {
1316  G4double PdotN = OldMomentum * theFacetNormal;
1317  G4double magP= OldMomentum.mag();
1318  G4double magN= theFacetNormal.mag();
1319  G4double incidentangle = pi - std::acos(PdotN/(magP*magN));
1320 
1321  return incidentangle;
1322 }
1323 
1324 G4double G4OpBoundaryProcess::GetReflectivity(G4double E1_perp,
1325  G4double E1_parl,
1326  G4double incidentangle,
1327  G4double RealRindex,
1328  G4double ImaginaryRindex)
1329 {
1330  G4complex Reflectivity, Reflectivity_TE, Reflectivity_TM;
1331  G4complex N1(Rindex1, 0), N2(RealRindex, ImaginaryRindex);
1332  G4complex CosPhi;
1333 
1334  G4complex u(1,0); //unit number 1
1335 
1336  G4complex numeratorTE; // E1_perp=1 E1_parl=0 -> TE polarization
1337  G4complex numeratorTM; // E1_parl=1 E1_perp=0 -> TM polarization
1338  G4complex denominatorTE, denominatorTM;
1339  G4complex rTM, rTE;
1340 
1341  G4MaterialPropertiesTable* aMaterialPropertiesTable =
1342  Material1->GetMaterialPropertiesTable();
1343  G4MaterialPropertyVector* aPropertyPointerR =
1344  aMaterialPropertiesTable->GetProperty(kREALRINDEX);
1345  G4MaterialPropertyVector* aPropertyPointerI =
1346  aMaterialPropertiesTable->GetProperty(kIMAGINARYRINDEX);
1347  if (aPropertyPointerR && aPropertyPointerI) {
1348  G4double RRindex = aPropertyPointerR->Value(thePhotonMomentum);
1349  G4double IRindex = aPropertyPointerI->Value(thePhotonMomentum);
1350  N1 = G4complex(RRindex,IRindex);
1351  }
1352 
1353  // Following two equations, rTM and rTE, are from: "Introduction To Modern
1354  // Optics" written by Fowles
1355 
1356  CosPhi=std::sqrt(u-((std::sin(incidentangle)*std::sin(incidentangle))*(N1*N1)/(N2*N2)));
1357 
1358  numeratorTE = N1*std::cos(incidentangle) - N2*CosPhi;
1359  denominatorTE = N1*std::cos(incidentangle) + N2*CosPhi;
1360  rTE = numeratorTE/denominatorTE;
1361 
1362  numeratorTM = N2*std::cos(incidentangle) - N1*CosPhi;
1363  denominatorTM = N2*std::cos(incidentangle) + N1*CosPhi;
1364  rTM = numeratorTM/denominatorTM;
1365 
1366  // This is my calculaton for reflectivity on a metalic surface
1367  // depending on the fraction of TE and TM polarization
1368  // when TE polarization, E1_parl=0 and E1_perp=1, R=abs(rTE)^2 and
1369  // when TM polarization, E1_parl=1 and E1_perp=0, R=abs(rTM)^2
1370 
1371  Reflectivity_TE = (rTE*conj(rTE))*(E1_perp*E1_perp)
1372  / (E1_perp*E1_perp + E1_parl*E1_parl);
1373  Reflectivity_TM = (rTM*conj(rTM))*(E1_parl*E1_parl)
1374  / (E1_perp*E1_perp + E1_parl*E1_parl);
1375  Reflectivity = Reflectivity_TE + Reflectivity_TM;
1376 
1377  do {
1378  if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TE))
1379  {iTE = -1;}else{iTE = 1;}
1380  if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TM))
1381  {iTM = -1;}else{iTM = 1;}
1382  // Loop checking, 13-Aug-2015, Peter Gumplinger
1383  } while(iTE<0&&iTM<0);
1384 
1385  return real(Reflectivity);
1386 
1387 }
1388 
1389 void G4OpBoundaryProcess::CalculateReflectivity()
1390 {
1391  G4double RealRindex =
1392  PropertyPointer1->Value(thePhotonMomentum);
1393  G4double ImaginaryRindex =
1394  PropertyPointer2->Value(thePhotonMomentum);
1395 
1396  // calculate FacetNormal
1397  if ( theFinish == ground ) {
1398  theFacetNormal =
1399  GetFacetNormal(OldMomentum, theGlobalNormal);
1400  } else {
1401  theFacetNormal = theGlobalNormal;
1402  }
1403 
1404  G4double PdotN = OldMomentum * theFacetNormal;
1405  cost1 = -PdotN;
1406 
1407  if (std::abs(cost1) < 1.0 - kCarTolerance) {
1408  sint1 = std::sqrt(1. - cost1*cost1);
1409  } else {
1410  sint1 = 0.0;
1411  }
1412 
1413  G4ThreeVector A_trans, A_paral, E1pp, E1pl;
1414  G4double E1_perp, E1_parl;
1415 
1416  if (sint1 > 0.0 ) {
1417  A_trans = OldMomentum.cross(theFacetNormal);
1418  A_trans = A_trans.unit();
1419  E1_perp = OldPolarization * A_trans;
1420  E1pp = E1_perp * A_trans;
1421  E1pl = OldPolarization - E1pp;
1422  E1_parl = E1pl.mag();
1423  }
1424  else {
1425  A_trans = OldPolarization;
1426  // Here we Follow Jackson's conventions and we set the
1427  // parallel component = 1 in case of a ray perpendicular
1428  // to the surface
1429  E1_perp = 0.0;
1430  E1_parl = 1.0;
1431  }
1432 
1433  //calculate incident angle
1434  G4double incidentangle = GetIncidentAngle();
1435 
1436  //calculate the reflectivity depending on incident angle,
1437  //polarization and complex refractive
1438 
1439  theReflectivity =
1440  GetReflectivity(E1_perp, E1_parl, incidentangle,
1441  RealRindex, ImaginaryRindex);
1442 }
1443 
1444 G4bool G4OpBoundaryProcess::InvokeSD(const G4Step* pStep)
1445 {
1446  G4Step aStep = *pStep;
1447 
1448  aStep.AddTotalEnergyDeposit(thePhotonMomentum);
1449 
1450  G4VSensitiveDetector* sd = aStep.GetPostStepPoint()->GetSensitiveDetector();
1451  if (sd) return sd->Hit(&aStep);
1452  else return false;
1453 }
G4VSensitiveDetector::Hit
G4bool Hit(G4Step *aStep)
Definition: G4VSensitiveDetector.hh:114
G4OpBoundaryProcess::theModel
G4OpticalSurfaceModel theModel
Definition: G4OpBoundaryProcess.hh:248
G4SwapObj
void G4SwapObj(T *a, T *b)
Definition: templates.hh:129
G4OpBoundaryProcessStatus
G4OpBoundaryProcessStatus
Definition: G4OpBoundaryProcess.hh:100
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4ParticleChange::ProposePolarization
void ProposePolarization(G4double Px, G4double Py, G4double Pz)
G4ExceptionDescription
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:45
G4OpticalSurface::GetPhiIndexMax
G4int GetPhiIndexMax(void) const
Definition: G4OpticalSurface.cc:228
G4VPhysicalVolume::GetLogicalVolume
G4LogicalVolume * GetLogicalVolume() const
G4OpticalSurface::GetThetaIndexMax
G4int GetThetaIndexMax(void) const
Definition: G4OpticalSurface.cc:223
G4Step::GetPreStepPoint
G4StepPoint * GetPreStepPoint() const
nm
static constexpr double nm
Definition: G4SIunits.hh:112
G4OpticalSurface::GetPolish
G4double GetPolish() const
Definition: G4OpticalSurface.hh:186
G4Track::GetStepLength
G4double GetStepLength() const
G4endl
#define G4endl
Definition: G4ios.hh:61
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
perCent
static constexpr double perCent
Definition: G4SIunits.hh:332
G4OpBoundaryProcess::GetReflectivity
G4double GetReflectivity(G4double E1_perp, G4double E1_parl, G4double incidentangle, G4double RealRindex, G4double ImaginaryRindex)
Definition: G4OpBoundaryProcess.cc:1324
G4StepPoint::GetSensitiveDetector
G4VSensitiveDetector * GetSensitiveDetector() const
G4ParticleChange::ProposeVelocity
void ProposeVelocity(G4double finalVelocity)
G4LogicalBorderSurface::GetSurface
static G4LogicalBorderSurface * GetSurface(const G4VPhysicalVolume *vol1, const G4VPhysicalVolume *vol2)
Definition: G4LogicalBorderSurface.cc:141
CLHEP::Hep3Vector::rotateUz
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
G4StepPoint::GetPhysicalVolume
G4VPhysicalVolume * GetPhysicalVolume() const
G4VDiscreteProcess::PostStepDoIt
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
Definition: G4VDiscreteProcess.cc:112
G4PhysicsVector::Value
G4double Value(G4double theEnergy, size_t &lastidx) const
Definition: G4PhysicsVector.cc:498
condition
G4double condition(const G4ErrorSymMatrix &m)
G4SurfaceType
G4SurfaceType
Definition: G4SurfaceProperty.hh:65
CLHEP::h_Planck
static constexpr double h_Planck
Definition: PhysicalConstants.h:63
G4OpBoundaryProcess::G4OpBoundaryProcess
G4OpBoundaryProcess(const G4String &processName="OpBoundary", G4ProcessType type=fOptical)
Definition: G4OpBoundaryProcess.cc:105
CLHEP::Hep3Vector::setX
void setX(double)
G4OpBoundaryProcess::PropertyPointer2
G4MaterialPropertyVector * PropertyPointer2
Definition: G4OpBoundaryProcess.hh:239
G4ParallelWorldProcess::GetHyperStep
static const G4Step * GetHyperStep()
Definition: G4ParallelWorldProcess.cc:56
G4OpBoundaryProcess::GetMeanFreePath
G4double GetMeanFreePath(const G4Track &, G4double, G4ForceCondition *condition)
Definition: G4OpBoundaryProcess.cc:1305
G4LogicalSkinSurface::GetSurface
static G4LogicalSkinSurface * GetSurface(const G4LogicalVolume *vol)
Definition: G4LogicalSkinSurface.cc:135
G4VProcess::aParticleChange
G4ParticleChange aParticleChange
Definition: G4VProcess.hh:289
CLHEP::Hep3Vector::angle
double angle(const Hep3Vector &) const
G4OpBoundaryProcess::DichroicVector
G4Physics2DVector * DichroicVector
Definition: G4OpBoundaryProcess.hh:265
G4complex
std::complex< G4double > G4complex
Definition: G4Types.hh:81
G4OpBoundaryProcess::PropertyPointer
G4MaterialPropertyVector * PropertyPointer
Definition: G4OpBoundaryProcess.hh:237
G4double
double G4double
Definition: G4Types.hh:76
G4bool
bool G4bool
Definition: G4Types.hh:79
G4StepPoint::GetStepStatus
G4StepStatus GetStepStatus() const
G4VPhysicalVolume::GetMotherLogical
G4LogicalVolume * GetMotherLogical() const
G4StepPoint::GetMaterial
G4Material * GetMaterial() const
G4OpticalSurface::GetDichroicVector
G4Physics2DVector * GetDichroicVector()
Definition: G4OpticalSurface.hh:295
G4ParticleChange::Initialize
virtual void Initialize(const G4Track &)
Definition: G4ParticleChange.cc:228
G4OpBoundaryProcess::theFinish
G4OpticalSurfaceFinish theFinish
Definition: G4OpBoundaryProcess.hh:250
G4OpticalSurface::GetModel
G4OpticalSurfaceModel GetModel() const
Definition: G4OpticalSurface.hh:175
CLHEP::Hep3Vector::setZ
void setZ(double)
G4ProcessType
G4ProcessType
Definition: G4ProcessType.hh:43
G4StepPoint::GetPosition
const G4ThreeVector & GetPosition() const
G4OpBoundaryProcess::OpticalSurface
G4OpticalSurface * OpticalSurface
Definition: G4OpBoundaryProcess.hh:235
alpha
static const G4double alpha
Definition: G4LivermoreBremsstrahlungModel.cc:79
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
twopi
static constexpr double twopi
Definition: G4SIunits.hh:76
G4VProcess::GetProcessName
const G4String & GetProcessName() const
Definition: G4VProcess.hh:411
halfpi
static constexpr double halfpi
Definition: G4SIunits.hh:77
G4LogicalSurface::GetSurfaceProperty
G4SurfaceProperty * GetSurfaceProperty() const
G4Step
Definition: G4Step.hh:76
G4Step::GetPostStepPoint
G4StepPoint * GetPostStepPoint() const
G4OpBoundaryProcess::theStatus
G4OpBoundaryProcessStatus theStatus
Definition: G4OpBoundaryProcess.hh:246
G4Track::GetVelocity
G4double GetVelocity() const
G4INCL::DeJongSpin::shoot
ThreeVector shoot(const G4int Ap, const G4int Af)
Definition: G4INCLDeJongSpin.cc:68
G4cerr
G4GLOB_DLL std::ostream G4cerr
CLHEP::Hep3Vector::rotate
Hep3Vector & rotate(double, const Hep3Vector &)
Definition: ThreeVectorR.cc:28
G4Step::AddTotalEnergyDeposit
void AddTotalEnergyDeposit(G4double value)
G4Physics2DVector::Value
G4double Value(G4double x, G4double y, size_t &lastidx, size_t &lastidy) const
Definition: G4Physics2DVector.cc:147
CLHEP::Hep3Vector::cross
Hep3Vector cross(const Hep3Vector &) const
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
G4TransportationManager::GetTransportationManager
static G4TransportationManager * GetTransportationManager()
Definition: G4TransportationManager.cc:96
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
G4OpticalSurface::GetSigmaAlpha
G4double GetSigmaAlpha() const
Definition: G4OpticalSurface.hh:181
CLHEP::Hep3Vector::mag
double mag() const
G4DynamicParticle::GetPolarization
const G4ThreeVector & GetPolarization() const
G4SwapPtr
void G4SwapPtr(T *&a, T *&b)
Definition: templates.hh:121
G4int
int G4int
Definition: G4Types.hh:78
CLHEP::c_light
static constexpr double c_light
Definition: PhysicalConstants.h:55
G4VProcess::verboseLevel
G4int verboseLevel
Definition: G4VProcess.hh:371
G4ParticleChange::ProposeMomentumDirection
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4OpBoundaryProcess::PostStepDoIt
G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep)
Definition: G4OpBoundaryProcess.cc:170
G4VProcess::SetProcessSubType
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:435
G4ForceCondition
G4ForceCondition
Definition: G4ForceCondition.hh:49
G4OpBoundaryProcess::G4BooleanRand
G4bool G4BooleanRand(const G4double prob) const
Definition: G4OpBoundaryProcess.hh:275
G4cout
G4GLOB_DLL std::ostream G4cout
G4MaterialPropertiesTable::GetProperty
G4MaterialPropertyVector * GetProperty(const char *key, G4bool warning=false)
Definition: G4MaterialPropertiesTable.cc:224
G4OpticalSurface::GetFinish
G4OpticalSurfaceFinish GetFinish() const
Definition: G4OpticalSurface.hh:170
G4GeometryTolerance::GetInstance
static G4GeometryTolerance * GetInstance()
Definition: G4GeometryTolerance.cc:73
CLHEP::Hep3Vector::orthogonal
Hep3Vector orthogonal() const
G4OpBoundaryProcess::PropertyPointer1
G4MaterialPropertyVector * PropertyPointer1
Definition: G4OpBoundaryProcess.hh:238
n
Char_t n[5]
Definition: comparison_ascii.C:55
G4DynamicParticle::GetTotalMomentum
G4double GetTotalMomentum() const
pi
static constexpr double pi
Definition: G4SIunits.hh:75
G4VParticleChange::ProposeLocalEnergyDeposit
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4SurfaceProperty::GetType
const G4SurfaceType & GetType() const
Definition: G4SurfaceProperty.hh:119
G4OpBoundaryProcess::GetFacetNormal
G4ThreeVector GetFacetNormal(const G4ThreeVector &Momentum, const G4ThreeVector &Normal) const
Definition: G4OpBoundaryProcess.cc:639
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83
G4VParticleChange::ProposeTrackStatus
void ProposeTrackStatus(G4TrackStatus status)
G4OpBoundaryProcess::InvokeSD
G4bool InvokeSD(const G4Step *step)
Definition: G4OpBoundaryProcess.cc:1444
G4Track::GetDynamicParticle
const G4DynamicParticle * GetDynamicParticle() const
G4GeometryTolerance::GetSurfaceTolerance
G4double GetSurfaceTolerance() const
Definition: G4GeometryTolerance.cc:87
CLHEP::Hep3Vector::setY
void setY(double)
G4OpBoundaryProcess::BoundaryProcessVerbose
void BoundaryProcessVerbose(void) const
Definition: G4OpBoundaryProcess.cc:554
G4VPhysicalVolume::GetName
const G4String & GetName() const
G4Material::GetMaterialPropertiesTable
G4MaterialPropertiesTable * GetMaterialPropertiesTable() const
Definition: G4Material.hh:252
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:117
G4OpticalSurface::GetLUTbins
G4int GetLUTbins(void) const
Definition: G4OpticalSurface.cc:213
다음에 의해 생성됨 :   doxygen 1.8.5