Geant4_10
G4OpBoundaryProcess.cc
Go to the documentation of this file.
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 //
69 // Author: Peter Gumplinger
70 // adopted from work by Werner Keil - April 2/96
71 // mail: gum@triumf.ca
72 //
74 
75 #include "G4ios.hh"
76 #include "G4PhysicalConstants.hh"
77 #include "G4OpProcessSubType.hh"
78 
79 #include "G4OpBoundaryProcess.hh"
80 #include "G4GeometryTolerance.hh"
81 
82 #include "G4VSensitiveDetector.hh"
84 
85 #include "G4SystemOfUnits.hh"
86 
88 // Class Implementation
90 
92  // Operators
94 
95 // G4OpBoundaryProcess::operator=(const G4OpBoundaryProcess &right)
96 // {
97 // }
98 
100  // Constructors
102 
104  G4ProcessType type)
105  : G4VDiscreteProcess(processName, type)
106 {
107  if ( verboseLevel > 0) {
108  G4cout << GetProcessName() << " is created " << G4endl;
109  }
110 
112 
113  theStatus = Undefined;
114  theModel = glisur;
115  theFinish = polished;
116  theReflectivity = 1.;
117  theEfficiency = 0.;
118  theTransmittance = 0.;
119 
120  prob_sl = 0.;
121  prob_ss = 0.;
122  prob_bs = 0.;
123 
124  PropertyPointer = NULL;
125  PropertyPointer1 = NULL;
126  PropertyPointer2 = NULL;
127 
128  Material1 = NULL;
129  Material2 = NULL;
130 
131  OpticalSurface = NULL;
132 
133  kCarTolerance = G4GeometryTolerance::GetInstance()
135 
136  iTE = iTM = 0;
137  thePhotonMomentum = 0.;
138  Rindex1 = Rindex2 = cost1 = cost2 = sint1 = sint2 = 0.;
139 
140  idx = idy = 0;
141  DichroicVector = NULL;
142 }
143 
144 // G4OpBoundaryProcess::G4OpBoundaryProcess(const G4OpBoundaryProcess &right)
145 // {
146 // }
147 
149  // Destructors
151 
153 
155  // Methods
157 
158 // PostStepDoIt
159 // ------------
160 //
161 
164 {
165  theStatus = Undefined;
166 
167  aParticleChange.Initialize(aTrack);
169 
170  // Get hyperStep from G4ParallelWorldProcess
171  // NOTE: PostSetpDoIt of this process should be
172  // invoked after G4ParallelWorldProcess!
173 
174  const G4Step* pStep = &aStep;
175 
177 
178  if (hStep) pStep = hStep;
179 
180  G4bool isOnBoundary =
182 
183  if (isOnBoundary) {
184  Material1 = pStep->GetPreStepPoint()->GetMaterial();
185  Material2 = pStep->GetPostStepPoint()->GetMaterial();
186  } else {
187  theStatus = NotAtBoundary;
188  if ( verboseLevel > 0) BoundaryProcessVerbose();
189  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
190  }
191 
192  G4VPhysicalVolume* thePrePV =
193  pStep->GetPreStepPoint() ->GetPhysicalVolume();
194  G4VPhysicalVolume* thePostPV =
196 
197  if ( verboseLevel > 0 ) {
198  G4cout << " Photon at Boundary! " << G4endl;
199  if (thePrePV) G4cout << " thePrePV: " << thePrePV->GetName() << G4endl;
200  if (thePostPV) G4cout << " thePostPV: " << thePostPV->GetName() << G4endl;
201  }
202 
203  if (aTrack.GetStepLength()<=kCarTolerance/2){
204  theStatus = StepTooSmall;
205  if ( verboseLevel > 0) BoundaryProcessVerbose();
206  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
207  }
208 
209  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
210 
211  thePhotonMomentum = aParticle->GetTotalMomentum();
212  OldMomentum = aParticle->GetMomentumDirection();
213  OldPolarization = aParticle->GetPolarization();
214 
215  if ( verboseLevel > 0 ) {
216  G4cout << " Old Momentum Direction: " << OldMomentum << G4endl;
217  G4cout << " Old Polarization: " << OldPolarization << G4endl;
218  }
219 
220  G4ThreeVector theGlobalPoint = pStep->GetPostStepPoint()->GetPosition();
221 
222  G4bool valid;
223  // Use the new method for Exit Normal in global coordinates,
224  // which provides the normal more reliably.
225 
226  // ID of Navigator which limits step
227 
229  std::vector<G4Navigator*>::iterator iNav =
231  GetActiveNavigatorsIterator();
232  theGlobalNormal =
233  (iNav[hNavId])->GetGlobalExitNormal(theGlobalPoint,&valid);
234 
235  if (valid) {
236  theGlobalNormal = -theGlobalNormal;
237  }
238  else
239  {
241  ed << " G4OpBoundaryProcess/PostStepDoIt(): "
242  << " The Navigator reports that it returned an invalid normal"
243  << G4endl;
244  G4Exception("G4OpBoundaryProcess::PostStepDoIt", "OpBoun01",
246  "Invalid Surface Normal - Geometry must return valid surface normal");
247  }
248 
249  if (OldMomentum * theGlobalNormal > 0.0) {
250 #ifdef G4OPTICAL_DEBUG
252  ed << " G4OpBoundaryProcess/PostStepDoIt(): "
253  << " theGlobalNormal points in a wrong direction. "
254  << G4endl;
255  ed << " The momentum of the photon arriving at interface (oldMomentum)"
256  << " must exit the volume cross in the step. " << G4endl;
257  ed << " So it MUST have dot < 0 with the normal that Exits the new volume (globalNormal)." << G4endl;
258  ed << " >> The dot product of oldMomentum and global Normal is " << OldMomentum*theGlobalNormal << G4endl;
259  ed << " Old Momentum (during step) = " << OldMomentum << G4endl;
260  ed << " Global Normal (Exiting New Vol) = " << theGlobalNormal << G4endl;
261  ed << G4endl;
262  G4Exception("G4OpBoundaryProcess::PostStepDoIt", "OpBoun02",
263  EventMustBeAborted, // Or JustWarning to see if it happens repeatedbly on one ray
264  ed,
265  "Invalid Surface Normal - Geometry must return valid surface normal pointing in the right direction");
266 #else
267  theGlobalNormal = -theGlobalNormal;
268 #endif
269  }
270 
271  G4MaterialPropertiesTable* aMaterialPropertiesTable;
272  G4MaterialPropertyVector* Rindex;
273 
274  aMaterialPropertiesTable = Material1->GetMaterialPropertiesTable();
275  if (aMaterialPropertiesTable) {
276  Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
277  }
278  else {
279  theStatus = NoRINDEX;
280  if ( verboseLevel > 0) BoundaryProcessVerbose();
281  aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
283  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
284  }
285 
286  if (Rindex) {
287  Rindex1 = Rindex->Value(thePhotonMomentum);
288  }
289  else {
290  theStatus = NoRINDEX;
291  if ( verboseLevel > 0) BoundaryProcessVerbose();
292  aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
294  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
295  }
296 
297  theReflectivity = 1.;
298  theEfficiency = 0.;
299  theTransmittance = 0.;
300 
301  theModel = glisur;
302  theFinish = polished;
303 
305 
306  Rindex = NULL;
307  OpticalSurface = NULL;
308 
309  G4LogicalSurface* Surface = NULL;
310 
311  Surface = G4LogicalBorderSurface::GetSurface(thePrePV, thePostPV);
312 
313  if (Surface == NULL){
314  G4bool enteredDaughter= (thePostPV->GetMotherLogical() ==
315  thePrePV ->GetLogicalVolume());
316  if(enteredDaughter){
317  Surface =
319  if(Surface == NULL)
320  Surface =
322  }
323  else {
324  Surface =
326  if(Surface == NULL)
327  Surface =
329  }
330  }
331 
332  if (Surface) OpticalSurface =
333  dynamic_cast <G4OpticalSurface*> (Surface->GetSurfaceProperty());
334 
335  if (OpticalSurface) {
336 
337  type = OpticalSurface->GetType();
338  theModel = OpticalSurface->GetModel();
339  theFinish = OpticalSurface->GetFinish();
340 
341  aMaterialPropertiesTable = OpticalSurface->
342  GetMaterialPropertiesTable();
343 
344  if (aMaterialPropertiesTable) {
345 
346  if (theFinish == polishedbackpainted ||
347  theFinish == groundbackpainted ) {
348  Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
349  if (Rindex) {
350  Rindex2 = Rindex->Value(thePhotonMomentum);
351  }
352  else {
353  theStatus = NoRINDEX;
354  if ( verboseLevel > 0) BoundaryProcessVerbose();
355  aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
357  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
358  }
359  }
360 
361  PropertyPointer =
362  aMaterialPropertiesTable->GetProperty("REFLECTIVITY");
363  PropertyPointer1 =
364  aMaterialPropertiesTable->GetProperty("REALRINDEX");
365  PropertyPointer2 =
366  aMaterialPropertiesTable->GetProperty("IMAGINARYRINDEX");
367 
368  iTE = 1;
369  iTM = 1;
370 
371  if (PropertyPointer) {
372 
373  theReflectivity =
374  PropertyPointer->Value(thePhotonMomentum);
375 
376  } else if (PropertyPointer1 && PropertyPointer2) {
377 
378  CalculateReflectivity();
379 
380  }
381 
382  PropertyPointer =
383  aMaterialPropertiesTable->GetProperty("EFFICIENCY");
384  if (PropertyPointer) {
385  theEfficiency =
386  PropertyPointer->Value(thePhotonMomentum);
387  }
388 
389  PropertyPointer =
390  aMaterialPropertiesTable->GetProperty("TRANSMITTANCE");
391  if (PropertyPointer) {
392  theTransmittance =
393  PropertyPointer->Value(thePhotonMomentum);
394  }
395 
396  if ( theModel == unified ) {
397  PropertyPointer =
398  aMaterialPropertiesTable->GetProperty("SPECULARLOBECONSTANT");
399  if (PropertyPointer) {
400  prob_sl =
401  PropertyPointer->Value(thePhotonMomentum);
402  } else {
403  prob_sl = 0.0;
404  }
405 
406  PropertyPointer =
407  aMaterialPropertiesTable->GetProperty("SPECULARSPIKECONSTANT");
408  if (PropertyPointer) {
409  prob_ss =
410  PropertyPointer->Value(thePhotonMomentum);
411  } else {
412  prob_ss = 0.0;
413  }
414 
415  PropertyPointer =
416  aMaterialPropertiesTable->GetProperty("BACKSCATTERCONSTANT");
417  if (PropertyPointer) {
418  prob_bs =
419  PropertyPointer->Value(thePhotonMomentum);
420  } else {
421  prob_bs = 0.0;
422  }
423  }
424  }
425  else if (theFinish == polishedbackpainted ||
426  theFinish == groundbackpainted ) {
427  aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
429  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
430  }
431  }
432 
433  if (type == dielectric_dielectric ) {
434  if (theFinish == polished || theFinish == ground ) {
435 
436  if (Material1 == Material2){
437  theStatus = SameMaterial;
438  if ( verboseLevel > 0) BoundaryProcessVerbose();
439  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
440  }
441  aMaterialPropertiesTable =
442  Material2->GetMaterialPropertiesTable();
443  if (aMaterialPropertiesTable)
444  Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
445  if (Rindex) {
446  Rindex2 = Rindex->Value(thePhotonMomentum);
447  }
448  else {
449  theStatus = NoRINDEX;
450  if ( verboseLevel > 0) BoundaryProcessVerbose();
451  aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
453  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
454  }
455  }
456  }
457 
458  if (type == dielectric_metal) {
459 
460  DielectricMetal();
461 
462  // Uncomment the following lines if you wish to have
463  // Transmission instead of Absorption
464  // if (theStatus == Absorption) {
465  // return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
466  // }
467 
468  }
469  else if (type == dielectric_LUT) {
470 
471  DielectricLUT();
472 
473  }
474  else if (type == dielectric_dichroic) {
475 
476  DielectricDichroic();
477 
478  }
479  else if (type == dielectric_dielectric) {
480 
481  if ( theFinish == polishedbackpainted ||
482  theFinish == groundbackpainted ) {
483  DielectricDielectric();
484  }
485  else {
486  if ( !G4BooleanRand(theReflectivity) ) {
487  DoAbsorption();
488  }
489  else {
490  if ( theFinish == polishedfrontpainted ) {
491  DoReflection();
492  }
493  else if ( theFinish == groundfrontpainted ) {
494  theStatus = LambertianReflection;
495  DoReflection();
496  }
497  else {
498  DielectricDielectric();
499  }
500  }
501  }
502  }
503  else {
504 
505  G4cerr << " Error: G4BoundaryProcess: illegal boundary type " << G4endl;
506  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
507 
508  }
509 
510  NewMomentum = NewMomentum.unit();
511  NewPolarization = NewPolarization.unit();
512 
513  if ( verboseLevel > 0) {
514  G4cout << " New Momentum Direction: " << NewMomentum << G4endl;
515  G4cout << " New Polarization: " << NewPolarization << G4endl;
516  BoundaryProcessVerbose();
517  }
518 
520  aParticleChange.ProposePolarization(NewPolarization);
521 
522  if ( theStatus == FresnelRefraction ) {
523  G4MaterialPropertyVector* groupvel =
524  Material2->GetMaterialPropertiesTable()->GetProperty("GROUPVEL");
525  G4double finalVelocity = groupvel->Value(thePhotonMomentum);
526  aParticleChange.ProposeVelocity(finalVelocity);
527  }
528 
529  if ( theStatus == Detection ) InvokeSD(pStep);
530 
531  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
532 }
533 
534 void G4OpBoundaryProcess::BoundaryProcessVerbose() const
535 {
536  if ( theStatus == Undefined )
537  G4cout << " *** Undefined *** " << G4endl;
538  if ( theStatus == FresnelRefraction )
539  G4cout << " *** FresnelRefraction *** " << G4endl;
540  if ( theStatus == FresnelReflection )
541  G4cout << " *** FresnelReflection *** " << G4endl;
542  if ( theStatus == TotalInternalReflection )
543  G4cout << " *** TotalInternalReflection *** " << G4endl;
544  if ( theStatus == LambertianReflection )
545  G4cout << " *** LambertianReflection *** " << G4endl;
546  if ( theStatus == LobeReflection )
547  G4cout << " *** LobeReflection *** " << G4endl;
548  if ( theStatus == SpikeReflection )
549  G4cout << " *** SpikeReflection *** " << G4endl;
550  if ( theStatus == BackScattering )
551  G4cout << " *** BackScattering *** " << G4endl;
552  if ( theStatus == PolishedLumirrorAirReflection )
553  G4cout << " *** PolishedLumirrorAirReflection *** " << G4endl;
554  if ( theStatus == PolishedLumirrorGlueReflection )
555  G4cout << " *** PolishedLumirrorGlueReflection *** " << G4endl;
556  if ( theStatus == PolishedAirReflection )
557  G4cout << " *** PolishedAirReflection *** " << G4endl;
558  if ( theStatus == PolishedTeflonAirReflection )
559  G4cout << " *** PolishedTeflonAirReflection *** " << G4endl;
560  if ( theStatus == PolishedTiOAirReflection )
561  G4cout << " *** PolishedTiOAirReflection *** " << G4endl;
562  if ( theStatus == PolishedTyvekAirReflection )
563  G4cout << " *** PolishedTyvekAirReflection *** " << G4endl;
564  if ( theStatus == PolishedVM2000AirReflection )
565  G4cout << " *** PolishedVM2000AirReflection *** " << G4endl;
566  if ( theStatus == PolishedVM2000GlueReflection )
567  G4cout << " *** PolishedVM2000GlueReflection *** " << G4endl;
568  if ( theStatus == EtchedLumirrorAirReflection )
569  G4cout << " *** EtchedLumirrorAirReflection *** " << G4endl;
570  if ( theStatus == EtchedLumirrorGlueReflection )
571  G4cout << " *** EtchedLumirrorGlueReflection *** " << G4endl;
572  if ( theStatus == EtchedAirReflection )
573  G4cout << " *** EtchedAirReflection *** " << G4endl;
574  if ( theStatus == EtchedTeflonAirReflection )
575  G4cout << " *** EtchedTeflonAirReflection *** " << G4endl;
576  if ( theStatus == EtchedTiOAirReflection )
577  G4cout << " *** EtchedTiOAirReflection *** " << G4endl;
578  if ( theStatus == EtchedTyvekAirReflection )
579  G4cout << " *** EtchedTyvekAirReflection *** " << G4endl;
580  if ( theStatus == EtchedVM2000AirReflection )
581  G4cout << " *** EtchedVM2000AirReflection *** " << G4endl;
582  if ( theStatus == EtchedVM2000GlueReflection )
583  G4cout << " *** EtchedVM2000GlueReflection *** " << G4endl;
584  if ( theStatus == GroundLumirrorAirReflection )
585  G4cout << " *** GroundLumirrorAirReflection *** " << G4endl;
586  if ( theStatus == GroundLumirrorGlueReflection )
587  G4cout << " *** GroundLumirrorGlueReflection *** " << G4endl;
588  if ( theStatus == GroundAirReflection )
589  G4cout << " *** GroundAirReflection *** " << G4endl;
590  if ( theStatus == GroundTeflonAirReflection )
591  G4cout << " *** GroundTeflonAirReflection *** " << G4endl;
592  if ( theStatus == GroundTiOAirReflection )
593  G4cout << " *** GroundTiOAirReflection *** " << G4endl;
594  if ( theStatus == GroundTyvekAirReflection )
595  G4cout << " *** GroundTyvekAirReflection *** " << G4endl;
596  if ( theStatus == GroundVM2000AirReflection )
597  G4cout << " *** GroundVM2000AirReflection *** " << G4endl;
598  if ( theStatus == GroundVM2000GlueReflection )
599  G4cout << " *** GroundVM2000GlueReflection *** " << G4endl;
600  if ( theStatus == Absorption )
601  G4cout << " *** Absorption *** " << G4endl;
602  if ( theStatus == Detection )
603  G4cout << " *** Detection *** " << G4endl;
604  if ( theStatus == NotAtBoundary )
605  G4cout << " *** NotAtBoundary *** " << G4endl;
606  if ( theStatus == SameMaterial )
607  G4cout << " *** SameMaterial *** " << G4endl;
608  if ( theStatus == StepTooSmall )
609  G4cout << " *** StepTooSmall *** " << G4endl;
610  if ( theStatus == NoRINDEX )
611  G4cout << " *** NoRINDEX *** " << G4endl;
612  if ( theStatus == Dichroic )
613  G4cout << " *** Dichroic Transmission *** " << G4endl;
614 }
615 
617 G4OpBoundaryProcess::GetFacetNormal(const G4ThreeVector& Momentum,
618  const G4ThreeVector& Normal ) const
619 {
620  G4ThreeVector FacetNormal;
621 
622  if (theModel == unified || theModel == LUT) {
623 
624  /* This function code alpha to a random value taken from the
625  distribution p(alpha) = g(alpha; 0, sigma_alpha)*std::sin(alpha),
626  for alpha > 0 and alpha < 90, where g(alpha; 0, sigma_alpha)
627  is a gaussian distribution with mean 0 and standard deviation
628  sigma_alpha. */
629 
630  G4double alpha;
631 
632  G4double sigma_alpha = 0.0;
633  if (OpticalSurface) sigma_alpha = OpticalSurface->GetSigmaAlpha();
634 
635  if (sigma_alpha == 0.0) return FacetNormal = Normal;
636 
637  G4double f_max = std::min(1.0,4.*sigma_alpha);
638 
639  do {
640  do {
641  alpha = G4RandGauss::shoot(0.0,sigma_alpha);
642  } while (G4UniformRand()*f_max > std::sin(alpha) || alpha >= halfpi );
643 
644  G4double phi = G4UniformRand()*twopi;
645 
646  G4double SinAlpha = std::sin(alpha);
647  G4double CosAlpha = std::cos(alpha);
648  G4double SinPhi = std::sin(phi);
649  G4double CosPhi = std::cos(phi);
650 
651  G4double unit_x = SinAlpha * CosPhi;
652  G4double unit_y = SinAlpha * SinPhi;
653  G4double unit_z = CosAlpha;
654 
655  FacetNormal.setX(unit_x);
656  FacetNormal.setY(unit_y);
657  FacetNormal.setZ(unit_z);
658 
659  G4ThreeVector tmpNormal = Normal;
660 
661  FacetNormal.rotateUz(tmpNormal);
662  } while (Momentum * FacetNormal >= 0.0);
663  }
664  else {
665 
666  G4double polish = 1.0;
667  if (OpticalSurface) polish = OpticalSurface->GetPolish();
668 
669  if (polish < 1.0) {
670  do {
671  G4ThreeVector smear;
672  do {
673  smear.setX(2.*G4UniformRand()-1.0);
674  smear.setY(2.*G4UniformRand()-1.0);
675  smear.setZ(2.*G4UniformRand()-1.0);
676  } while (smear.mag()>1.0);
677  smear = (1.-polish) * smear;
678  FacetNormal = Normal + smear;
679  } while (Momentum * FacetNormal >= 0.0);
680  FacetNormal = FacetNormal.unit();
681  }
682  else {
683  FacetNormal = Normal;
684  }
685  }
686  return FacetNormal;
687 }
688 
689 void G4OpBoundaryProcess::DielectricMetal()
690 {
691  G4int n = 0;
692 
693  do {
694 
695  n++;
696 
697  if( !G4BooleanRand(theReflectivity) && n == 1 ) {
698 
699  // Comment out DoAbsorption and uncomment theStatus = Absorption;
700  // if you wish to have Transmission instead of Absorption
701 
702  DoAbsorption();
703  // theStatus = Absorption;
704  break;
705 
706  }
707  else {
708 
709  if (PropertyPointer1 && PropertyPointer2) {
710  if ( n > 1 ) {
711  CalculateReflectivity();
712  if ( !G4BooleanRand(theReflectivity) ) {
713  DoAbsorption();
714  break;
715  }
716  }
717  }
718 
719  if ( theModel == glisur || theFinish == polished ) {
720 
721  DoReflection();
722 
723  } else {
724 
725  if ( n == 1 ) ChooseReflection();
726 
727  if ( theStatus == LambertianReflection ) {
728  DoReflection();
729  }
730  else if ( theStatus == BackScattering ) {
731  NewMomentum = -OldMomentum;
732  NewPolarization = -OldPolarization;
733  }
734  else {
735 
736  if(theStatus==LobeReflection){
737  if ( PropertyPointer1 && PropertyPointer2 ){
738  } else {
739  theFacetNormal =
740  GetFacetNormal(OldMomentum,theGlobalNormal);
741  }
742  }
743 
744  G4double PdotN = OldMomentum * theFacetNormal;
745  NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
746  G4double EdotN = OldPolarization * theFacetNormal;
747 
748  G4ThreeVector A_trans, A_paral;
749 
750  if (sint1 > 0.0 ) {
751  A_trans = OldMomentum.cross(theFacetNormal);
752  A_trans = A_trans.unit();
753  } else {
754  A_trans = OldPolarization;
755  }
756  A_paral = NewMomentum.cross(A_trans);
757  A_paral = A_paral.unit();
758 
759  if(iTE>0&&iTM>0) {
760  NewPolarization =
761  -OldPolarization + (2.*EdotN)*theFacetNormal;
762  } else if (iTE>0) {
763  NewPolarization = -A_trans;
764  } else if (iTM>0) {
765  NewPolarization = -A_paral;
766  }
767 
768  }
769 
770  }
771 
772  OldMomentum = NewMomentum;
773  OldPolarization = NewPolarization;
774 
775  }
776 
777  } while (NewMomentum * theGlobalNormal < 0.0);
778 }
779 
780 void G4OpBoundaryProcess::DielectricLUT()
781 {
782  G4int thetaIndex, phiIndex;
783  G4double AngularDistributionValue, thetaRad, phiRad, EdotN;
784  G4ThreeVector PerpendicularVectorTheta, PerpendicularVectorPhi;
785 
786  theStatus = G4OpBoundaryProcessStatus(G4int(theFinish) +
788 
789  G4int thetaIndexMax = OpticalSurface->GetThetaIndexMax();
790  G4int phiIndexMax = OpticalSurface->GetPhiIndexMax();
791 
792  do {
793  if ( !G4BooleanRand(theReflectivity) ) // Not reflected, so Absorbed
794  DoAbsorption();
795  else {
796  // Calculate Angle between Normal and Photon Momentum
797  G4double anglePhotonToNormal =
798  OldMomentum.angle(-theGlobalNormal);
799  // Round it to closest integer
800  G4int angleIncident = G4int(std::floor(180/pi*anglePhotonToNormal+0.5));
801 
802  // Take random angles THETA and PHI,
803  // and see if below Probability - if not - Redo
804  do {
805  thetaIndex = G4RandFlat::shootInt(thetaIndexMax-1);
806  phiIndex = G4RandFlat::shootInt(phiIndexMax-1);
807  // Find probability with the new indeces from LUT
808  AngularDistributionValue = OpticalSurface ->
809  GetAngularDistributionValue(angleIncident,
810  thetaIndex,
811  phiIndex);
812  } while ( !G4BooleanRand(AngularDistributionValue) );
813 
814  thetaRad = (-90 + 4*thetaIndex)*pi/180;
815  phiRad = (-90 + 5*phiIndex)*pi/180;
816  // Rotate Photon Momentum in Theta, then in Phi
817  NewMomentum = -OldMomentum;
818  PerpendicularVectorTheta = NewMomentum.cross(theGlobalNormal);
819  if (PerpendicularVectorTheta.mag() > kCarTolerance ) {
820  PerpendicularVectorPhi =
821  PerpendicularVectorTheta.cross(NewMomentum);
822  }
823  else {
824  PerpendicularVectorTheta = NewMomentum.orthogonal();
825  PerpendicularVectorPhi =
826  PerpendicularVectorTheta.cross(NewMomentum);
827  }
828  NewMomentum =
829  NewMomentum.rotate(anglePhotonToNormal-thetaRad,
830  PerpendicularVectorTheta);
831  NewMomentum = NewMomentum.rotate(-phiRad,PerpendicularVectorPhi);
832  // Rotate Polarization too:
833  theFacetNormal = (NewMomentum - OldMomentum).unit();
834  EdotN = OldPolarization * theFacetNormal;
835  NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
836  }
837  } while (NewMomentum * theGlobalNormal <= 0.0);
838 }
839 
840 void G4OpBoundaryProcess::DielectricDichroic()
841 {
842  // Calculate Angle between Normal and Photon Momentum
843  G4double anglePhotonToNormal = OldMomentum.angle(-theGlobalNormal);
844 
845  // Round it to closest integer
846  G4double angleIncident = std::floor(180/pi*anglePhotonToNormal+0.5);
847 
848  if (!DichroicVector) {
849  if (OpticalSurface) DichroicVector = OpticalSurface->GetDichroicVector();
850  }
851 
852 
853  if (DichroicVector) {
854  G4double wavelength = h_Planck*c_light/thePhotonMomentum;
855  theTransmittance =
856  DichroicVector->Value(wavelength/nm,angleIncident,idx,idy)*perCent;
857 // G4cout << "wavelength: " << std::floor(wavelength/nm)
858 // << "nm" << G4endl;
859 // G4cout << "Incident angle: " << angleIncident << "deg" << G4endl;
860 // G4cout << "Transmittance: "
861 // << std::floor(theTransmittance/perCent) << "%" << G4endl;
862  } else {
864  ed << " G4OpBoundaryProcess/DielectricDichroic(): "
865  << " The dichroic surface has no G4Physics2DVector"
866  << G4endl;
867  G4Exception("G4OpBoundaryProcess::DielectricDichroic", "OpBoun03",
868  FatalException,ed,
869  "A dichroic surface must have an associated G4Physics2DVector");
870  }
871 
872  if ( !G4BooleanRand(theTransmittance) ) // Not transmitted, so reflect
873  DoReflection();
874  else {
875  theStatus = Dichroic;
876  NewMomentum = OldMomentum;
877  NewPolarization = OldPolarization;
878  }
879 }
880 
881 void G4OpBoundaryProcess::DielectricDielectric()
882 {
883  G4bool Inside = false;
884  G4bool Swap = false;
885 
886  leap:
887 
888  G4bool Through = false;
889  G4bool Done = false;
890 
891  do {
892 
893  if (Through) {
894  Swap = !Swap;
895  Through = false;
896  theGlobalNormal = -theGlobalNormal;
897  G4SwapPtr(Material1,Material2);
898  G4SwapObj(&Rindex1,&Rindex2);
899  }
900 
901  if ( theFinish == polished ) {
902  theFacetNormal = theGlobalNormal;
903  }
904  else {
905  theFacetNormal =
906  GetFacetNormal(OldMomentum,theGlobalNormal);
907  }
908 
909  G4double PdotN = OldMomentum * theFacetNormal;
910  G4double EdotN = OldPolarization * theFacetNormal;
911 
912  cost1 = - PdotN;
913  if (std::abs(cost1) < 1.0-kCarTolerance){
914  sint1 = std::sqrt(1.-cost1*cost1);
915  sint2 = sint1*Rindex1/Rindex2; // *** Snell's Law ***
916  }
917  else {
918  sint1 = 0.0;
919  sint2 = 0.0;
920  }
921 
922  if (sint2 >= 1.0) {
923 
924  // Simulate total internal reflection
925 
926  if (Swap) Swap = !Swap;
927 
928  theStatus = TotalInternalReflection;
929 
930  if ( theModel == unified && theFinish != polished )
931  ChooseReflection();
932 
933  if ( theStatus == LambertianReflection ) {
934  DoReflection();
935  }
936  else if ( theStatus == BackScattering ) {
937  NewMomentum = -OldMomentum;
938  NewPolarization = -OldPolarization;
939  }
940  else {
941 
942  PdotN = OldMomentum * theFacetNormal;
943  NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
944  EdotN = OldPolarization * theFacetNormal;
945  NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
946 
947  }
948  }
949  else if (sint2 < 1.0) {
950 
951  // Calculate amplitude for transmission (Q = P x N)
952 
953  if (cost1 > 0.0) {
954  cost2 = std::sqrt(1.-sint2*sint2);
955  }
956  else {
957  cost2 = -std::sqrt(1.-sint2*sint2);
958  }
959 
960  G4ThreeVector A_trans, A_paral, E1pp, E1pl;
961  G4double E1_perp, E1_parl;
962 
963  if (sint1 > 0.0) {
964  A_trans = OldMomentum.cross(theFacetNormal);
965  A_trans = A_trans.unit();
966  E1_perp = OldPolarization * A_trans;
967  E1pp = E1_perp * A_trans;
968  E1pl = OldPolarization - E1pp;
969  E1_parl = E1pl.mag();
970  }
971  else {
972  A_trans = OldPolarization;
973  // Here we Follow Jackson's conventions and we set the
974  // parallel component = 1 in case of a ray perpendicular
975  // to the surface
976  E1_perp = 0.0;
977  E1_parl = 1.0;
978  }
979 
980  G4double s1 = Rindex1*cost1;
981  G4double E2_perp = 2.*s1*E1_perp/(Rindex1*cost1+Rindex2*cost2);
982  G4double E2_parl = 2.*s1*E1_parl/(Rindex2*cost1+Rindex1*cost2);
983  G4double E2_total = E2_perp*E2_perp + E2_parl*E2_parl;
984  G4double s2 = Rindex2*cost2*E2_total;
985 
986  G4double TransCoeff;
987 
988  if (theTransmittance > 0) TransCoeff = theTransmittance;
989  else if (cost1 != 0.0) TransCoeff = s2/s1;
990  else TransCoeff = 0.0;
991 
992  G4double E2_abs, C_parl, C_perp;
993 
994  if ( !G4BooleanRand(TransCoeff) ) {
995 
996  // Simulate reflection
997 
998  if (Swap) Swap = !Swap;
999 
1000  theStatus = FresnelReflection;
1001 
1002  if ( theModel == unified && theFinish != polished )
1003  ChooseReflection();
1004 
1005  if ( theStatus == LambertianReflection ) {
1006  DoReflection();
1007  }
1008  else if ( theStatus == BackScattering ) {
1009  NewMomentum = -OldMomentum;
1010  NewPolarization = -OldPolarization;
1011  }
1012  else {
1013 
1014  PdotN = OldMomentum * theFacetNormal;
1015  NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
1016 
1017  if (sint1 > 0.0) { // incident ray oblique
1018 
1019  E2_parl = Rindex2*E2_parl/Rindex1 - E1_parl;
1020  E2_perp = E2_perp - E1_perp;
1021  E2_total = E2_perp*E2_perp + E2_parl*E2_parl;
1022  A_paral = NewMomentum.cross(A_trans);
1023  A_paral = A_paral.unit();
1024  E2_abs = std::sqrt(E2_total);
1025  C_parl = E2_parl/E2_abs;
1026  C_perp = E2_perp/E2_abs;
1027 
1028  NewPolarization = C_parl*A_paral + C_perp*A_trans;
1029 
1030  }
1031 
1032  else { // incident ray perpendicular
1033 
1034  if (Rindex2 > Rindex1) {
1035  NewPolarization = - OldPolarization;
1036  }
1037  else {
1038  NewPolarization = OldPolarization;
1039  }
1040 
1041  }
1042  }
1043  }
1044  else { // photon gets transmitted
1045 
1046  // Simulate transmission/refraction
1047 
1048  Inside = !Inside;
1049  Through = true;
1050  theStatus = FresnelRefraction;
1051 
1052  if (sint1 > 0.0) { // incident ray oblique
1053 
1054  G4double alpha = cost1 - cost2*(Rindex2/Rindex1);
1055  NewMomentum = OldMomentum + alpha*theFacetNormal;
1056  NewMomentum = NewMomentum.unit();
1057  PdotN = -cost2;
1058  A_paral = NewMomentum.cross(A_trans);
1059  A_paral = A_paral.unit();
1060  E2_abs = std::sqrt(E2_total);
1061  C_parl = E2_parl/E2_abs;
1062  C_perp = E2_perp/E2_abs;
1063 
1064  NewPolarization = C_parl*A_paral + C_perp*A_trans;
1065 
1066  }
1067  else { // incident ray perpendicular
1068 
1069  NewMomentum = OldMomentum;
1070  NewPolarization = OldPolarization;
1071 
1072  }
1073  }
1074  }
1075 
1076  OldMomentum = NewMomentum.unit();
1077  OldPolarization = NewPolarization.unit();
1078 
1079  if (theStatus == FresnelRefraction) {
1080  Done = (NewMomentum * theGlobalNormal <= 0.0);
1081  }
1082  else {
1083  Done = (NewMomentum * theGlobalNormal >= 0.0);
1084  }
1085 
1086  } while (!Done);
1087 
1088  if (Inside && !Swap) {
1089  if( theFinish == polishedbackpainted ||
1090  theFinish == groundbackpainted ) {
1091 
1092  if( !G4BooleanRand(theReflectivity) ) {
1093  DoAbsorption();
1094  }
1095  else {
1096  if (theStatus != FresnelRefraction ) {
1097  theGlobalNormal = -theGlobalNormal;
1098  }
1099  else {
1100  Swap = !Swap;
1101  G4SwapPtr(Material1,Material2);
1102  G4SwapObj(&Rindex1,&Rindex2);
1103  }
1104  if ( theFinish == groundbackpainted )
1105  theStatus = LambertianReflection;
1106 
1107  DoReflection();
1108 
1109  theGlobalNormal = -theGlobalNormal;
1110  OldMomentum = NewMomentum;
1111 
1112  goto leap;
1113  }
1114  }
1115  }
1116 }
1117 
1118 // GetMeanFreePath
1119 // ---------------
1120 //
1122  G4double ,
1124 {
1125  *condition = Forced;
1126 
1127  return DBL_MAX;
1128 }
1129 
1130 G4double G4OpBoundaryProcess::GetIncidentAngle()
1131 {
1132  G4double PdotN = OldMomentum * theFacetNormal;
1133  G4double magP= OldMomentum.mag();
1134  G4double magN= theFacetNormal.mag();
1135  G4double incidentangle = pi - std::acos(PdotN/(magP*magN));
1136 
1137  return incidentangle;
1138 }
1139 
1140 G4double G4OpBoundaryProcess::GetReflectivity(G4double E1_perp,
1141  G4double E1_parl,
1142  G4double incidentangle,
1143  G4double RealRindex,
1144  G4double ImaginaryRindex)
1145 {
1146 
1147  G4complex Reflectivity, Reflectivity_TE, Reflectivity_TM;
1148  G4complex N(RealRindex, ImaginaryRindex);
1149  G4complex CosPhi;
1150 
1151  G4complex u(1,0); //unit number 1
1152 
1153  G4complex numeratorTE; // E1_perp=1 E1_parl=0 -> TE polarization
1154  G4complex numeratorTM; // E1_parl=1 E1_perp=0 -> TM polarization
1155  G4complex denominatorTE, denominatorTM;
1156  G4complex rTM, rTE;
1157 
1158  // Following two equations, rTM and rTE, are from: "Introduction To Modern
1159  // Optics" written by Fowles
1160 
1161  CosPhi=std::sqrt(u-((std::sin(incidentangle)*std::sin(incidentangle))/(N*N)));
1162 
1163  numeratorTE = std::cos(incidentangle) - N*CosPhi;
1164  denominatorTE = std::cos(incidentangle) + N*CosPhi;
1165  rTE = numeratorTE/denominatorTE;
1166 
1167  numeratorTM = N*std::cos(incidentangle) - CosPhi;
1168  denominatorTM = N*std::cos(incidentangle) + CosPhi;
1169  rTM = numeratorTM/denominatorTM;
1170 
1171  // This is my calculaton for reflectivity on a metalic surface
1172  // depending on the fraction of TE and TM polarization
1173  // when TE polarization, E1_parl=0 and E1_perp=1, R=abs(rTE)^2 and
1174  // when TM polarization, E1_parl=1 and E1_perp=0, R=abs(rTM)^2
1175 
1176  Reflectivity_TE = (rTE*conj(rTE))*(E1_perp*E1_perp)
1177  / (E1_perp*E1_perp + E1_parl*E1_parl);
1178  Reflectivity_TM = (rTM*conj(rTM))*(E1_parl*E1_parl)
1179  / (E1_perp*E1_perp + E1_parl*E1_parl);
1180  Reflectivity = Reflectivity_TE + Reflectivity_TM;
1181 
1182  do {
1183  if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TE))
1184  {iTE = -1;}else{iTE = 1;}
1185  if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TM))
1186  {iTM = -1;}else{iTM = 1;}
1187  } while(iTE<0&&iTM<0);
1188 
1189  return real(Reflectivity);
1190 
1191 }
1192 
1193 void G4OpBoundaryProcess::CalculateReflectivity()
1194 {
1195  G4double RealRindex =
1196  PropertyPointer1->Value(thePhotonMomentum);
1197  G4double ImaginaryRindex =
1198  PropertyPointer2->Value(thePhotonMomentum);
1199 
1200  // calculate FacetNormal
1201  if ( theFinish == ground ) {
1202  theFacetNormal =
1203  GetFacetNormal(OldMomentum, theGlobalNormal);
1204  } else {
1205  theFacetNormal = theGlobalNormal;
1206  }
1207 
1208  G4double PdotN = OldMomentum * theFacetNormal;
1209  cost1 = -PdotN;
1210 
1211  if (std::abs(cost1) < 1.0 - kCarTolerance) {
1212  sint1 = std::sqrt(1. - cost1*cost1);
1213  } else {
1214  sint1 = 0.0;
1215  }
1216 
1217  G4ThreeVector A_trans, A_paral, E1pp, E1pl;
1218  G4double E1_perp, E1_parl;
1219 
1220  if (sint1 > 0.0 ) {
1221  A_trans = OldMomentum.cross(theFacetNormal);
1222  A_trans = A_trans.unit();
1223  E1_perp = OldPolarization * A_trans;
1224  E1pp = E1_perp * A_trans;
1225  E1pl = OldPolarization - E1pp;
1226  E1_parl = E1pl.mag();
1227  }
1228  else {
1229  A_trans = OldPolarization;
1230  // Here we Follow Jackson's conventions and we set the
1231  // parallel component = 1 in case of a ray perpendicular
1232  // to the surface
1233  E1_perp = 0.0;
1234  E1_parl = 1.0;
1235  }
1236 
1237  //calculate incident angle
1238  G4double incidentangle = GetIncidentAngle();
1239 
1240  //calculate the reflectivity depending on incident angle,
1241  //polarization and complex refractive
1242 
1243  theReflectivity =
1244  GetReflectivity(E1_perp, E1_parl, incidentangle,
1245  RealRindex, ImaginaryRindex);
1246 }
1247 
1248 G4bool G4OpBoundaryProcess::InvokeSD(const G4Step* pStep)
1249 {
1250  G4Step aStep = *pStep;
1251 
1252  aStep.AddTotalEnergyDeposit(thePhotonMomentum);
1253 
1255  if (sd) return sd->Hit(&aStep);
1256  else return false;
1257 }
G4double condition(const G4ErrorSymMatrix &m)
ThreeVector shoot(const G4int Ap, const G4int Af)
void G4SwapObj(T *a, T *b)
Definition: templates.hh:129
G4MaterialPropertyVector * GetProperty(const char *key)
G4OpticalSurfaceModel GetModel() const
G4int verboseLevel
Definition: G4VProcess.hh:368
double angle(const Hep3Vector &) const
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double GetVelocity() const
const G4DynamicParticle * GetDynamicParticle() const
float h_Planck
Definition: hepunit.py:263
G4StepStatus GetStepStatus() const
G4Material * GetMaterial() const
G4double GetSurfaceTolerance() const
G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep)
void ProposePolarization(G4double Px, G4double Py, G4double Pz)
const G4SurfaceType & GetType() const
G4OpBoundaryProcessStatus
G4OpBoundaryProcess(const G4String &processName="OpBoundary", G4ProcessType type=fOptical)
G4SurfaceType
int G4int
Definition: G4Types.hh:78
void setY(double)
void setZ(double)
void setX(double)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
static G4LogicalBorderSurface * GetSurface(const G4VPhysicalVolume *vol1, const G4VPhysicalVolume *vol2)
G4double GetTotalMomentum() const
G4StepPoint * GetPreStepPoint() const
Char_t n[5]
std::complex< G4double > G4complex
Definition: G4Types.hh:81
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
G4double Value(G4double x, G4double y, size_t &lastidx, size_t &lastidy) const
G4double GetPolish() const
G4VPhysicalVolume * GetPhysicalVolume() const
const G4String & GetName() const
const G4ThreeVector & GetPosition() const
G4int GetPhiIndexMax(void) const
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:72
G4bool Hit(G4Step *aStep)
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:432
G4double GetMeanFreePath(const G4Track &, G4double, G4ForceCondition *condition)
Definition: G4Step.hh:76
G4LogicalVolume * GetMotherLogical() const
G4double Value(G4double theEnergy, size_t &lastidx) const
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
G4OpticalSurfaceFinish GetFinish() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static G4TransportationManager * GetTransportationManager()
virtual void Initialize(const G4Track &)
G4LogicalVolume * GetLogicalVolume() const
G4Physics2DVector * GetDichroicVector()
void G4SwapPtr(T *&a, T *&b)
Definition: templates.hh:121
G4MaterialPropertiesTable * GetMaterialPropertiesTable() const
Definition: G4Material.hh:250
Hep3Vector unit() const
G4StepPoint * GetPostStepPoint() const
void AddTotalEnergyDeposit(G4double value)
Hep3Vector orthogonal() const
const G4ThreeVector & GetPolarization() const
G4ParticleChange aParticleChange
Definition: G4VProcess.hh:289
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
float perCent
Definition: hepunit.py:239
#define G4endl
Definition: G4ios.hh:61
**D E S C R I P T I O N
Definition: HEPEvtcom.cc:77
G4SurfaceProperty * GetSurfaceProperty() const
Hep3Vector cross(const Hep3Vector &) const
G4VSensitiveDetector * GetSensitiveDetector() const
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
G4double GetSigmaAlpha() const
G4ForceCondition
void ProposeVelocity(G4double finalVelocity)
double mag() const
#define DBL_MAX
Definition: templates.hh:83
G4int GetThetaIndexMax(void) const
Hep3Vector & rotate(double, const Hep3Vector &)
Definition: ThreeVectorR.cc:28
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
float c_light
Definition: hepunit.py:257
G4double GetStepLength() const
static G4GeometryTolerance * GetInstance()
static const G4Step * GetHyperStep()
G4GLOB_DLL std::ostream G4cerr
G4ProcessType
static G4LogicalSkinSurface * GetSurface(const G4LogicalVolume *vol)