ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4OpMieHG.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4OpMieHG.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 //
26 //
28 //
29 // File G4OpMieHG.hh
30 // Description: Discrete Process -- Mie Scattering of Optical Photons
31 // Created: 2010-07-03
32 // Author: Xin Qian
33 // Based on work from Vlasios Vasileiou
34 //
35 // This subroutine will mimic the Mie scattering based on
36 // Henyey-Greenstein phase function
37 // Forward and backward angles are treated separately.
38 //
39 // mail: gum@triumf.ca
40 //
42 
43 #include "G4OpMieHG.hh"
44 #include "G4PhysicalConstants.hh"
45 #include "G4OpProcessSubType.hh"
46 
47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
48 
49 G4OpMieHG::G4OpMieHG(const G4String& processName, G4ProcessType type)
50  : G4VDiscreteProcess(processName, type)
51 {
52  if (verboseLevel>0) {
53  G4cout << GetProcessName() << " is created " << G4endl;
54  }
55 
57 }
58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
59 
61 
62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
63 
65 G4OpMieHG::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
66 {
68 
69  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
70  const G4Material* aMaterial = aTrack.GetMaterial();
71  G4MaterialPropertiesTable* aMaterialPropertyTable =
72  aMaterial->GetMaterialPropertiesTable();
73 
74  G4double forward_g =
75  aMaterialPropertyTable->GetConstProperty(kMIEHG_FORWARD);
76  G4double backward_g =
77  aMaterialPropertyTable->GetConstProperty(kMIEHG_BACKWARD);
78  G4double ForwardRatio =
79  aMaterialPropertyTable->GetConstProperty(kMIEHG_FORWARD_RATIO);
80 
81  if (verboseLevel >0 ) {
82  G4cout << "MIE Scattering Photon!" << G4endl;
83  G4cout << "MIE Old Momentum Direction: "
84  << aParticle->GetMomentumDirection() << G4endl;
85  G4cout << "MIE Old Polarization: "
86  << aParticle->GetPolarization() << G4endl;
87  }
88 
89  G4double gg;
90  G4int direction;
91  if (G4UniformRand() <= ForwardRatio){
92  gg = forward_g;
93  direction = 1;
94  } else {
95  gg = backward_g;
96  direction = -1;
97  }
98 
100 
101  G4double Theta;
102  //sample the direction
103  if (gg != 0.) {
104  Theta = std::acos(2.*r*(1.+gg)*(1.+gg)*(1.-gg+gg*r)/((1.-gg+2.*gg*r)*(1.-gg+2.*gg*r)) -1.);
105  } else {
106  Theta = std::acos(2.*r-1.);
107  }
108  G4double Phi = G4UniformRand()*twopi;
109  //G4double Phi = G4UniformRand()*2*pi;
110 
111  if (direction == -1) Theta = pi - Theta; //backward scattering
112 
113  G4ThreeVector NewMomentumDirection, OldMomentumDirection;
114  G4ThreeVector OldPolarization, NewPolarization;
115 
116  NewMomentumDirection.set
117  (std::sin(Theta)*std::cos(Phi), std::sin(Theta)*std::sin(Phi), std::cos(Theta));
118  OldMomentumDirection = aParticle->GetMomentumDirection();
119  NewMomentumDirection.rotateUz(OldMomentumDirection);
120  NewMomentumDirection = NewMomentumDirection.unit();
121 
122  OldPolarization = aParticle->GetPolarization();
123  G4double constant = -1./NewMomentumDirection.dot(OldPolarization);
124 
125  NewPolarization = NewMomentumDirection + constant*OldPolarization;
126  NewPolarization = NewPolarization.unit();
127 
128  if (NewPolarization.mag() == 0.) {
129  r = G4UniformRand()*twopi;
130  NewPolarization.set(std::cos(r),std::sin(r),0.);
131  NewPolarization.rotateUz(NewMomentumDirection);
132  } else {
133  // There are two directions which perpendicular
134  // new momentum direction
135  if (G4UniformRand() < 0.5) NewPolarization = -NewPolarization;
136  }
137 
138  aParticleChange.ProposePolarization(NewPolarization);
139  aParticleChange.ProposeMomentumDirection(NewMomentumDirection);
140 
141  if (verboseLevel > 0) {
142  G4cout << "MIE New Polarization: " << NewPolarization << G4endl;
143  G4cout << "MIE Polarization Change: " << *(aParticleChange.GetPolarization()) << G4endl;
144  G4cout << "MIE New Momentum Direction: " << NewMomentumDirection << G4endl;
145  G4cout << "MIE Momentum Change: " << *(aParticleChange.GetMomentumDirection()) << G4endl;
146  }
147 
148  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
149 }
150 
151 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
152 
154  G4double,
156 {
157  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
158  const G4Material* aMaterial = aTrack.GetMaterial();
159 
160  G4double thePhotonEnergy = aParticle->GetTotalEnergy();
161 
162  G4double AttenuationLength = DBL_MAX;
163 
164  G4MaterialPropertiesTable* aMaterialPropertyTable =
165  aMaterial->GetMaterialPropertiesTable();
166 
167  if (aMaterialPropertyTable) {
168  G4MaterialPropertyVector* AttenuationLengthVector =
169  aMaterialPropertyTable->GetProperty(kMIEHG);
170  if (AttenuationLengthVector) {
171  AttenuationLength = AttenuationLengthVector->Value(thePhotonEnergy);
172  }
173  // else {
174  // G4cout << "No Mie scattering length specified" << G4endl;
175  // }
176  }
177  //else {
178  // G4cout << "No Mie scattering length specified" << G4endl;
179  // }
180 
181 // G4cout << thePhotonEnergy/GeV << " \t" << AttenuationLength/m << G4endl;
182 
183  return AttenuationLength;
184 }