Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4OpWLS.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//
26//
27//
28////////////////////////////////////////////////////////////////////////
29// Optical Photon WaveLength Shifting (WLS) Class Implementation
30////////////////////////////////////////////////////////////////////////
31//
32// File: G4OpWLS.cc
33// Description: Discrete Process -- Wavelength Shifting of Optical Photons
34// Version: 1.0
35// Created: 2003-05-13
36// Author: John Paul Archambault
37// (Adaptation of G4Scintillation and G4OpAbsorption)
38// Updated: 2005-07-28 - add G4ProcessType to constructor
39// 2006-05-07 - add G4VWLSTimeGeneratorProfile
40//
41////////////////////////////////////////////////////////////////////////
42
43#include "G4OpWLS.hh"
44#include "G4ios.hh"
46#include "G4SystemOfUnits.hh"
47#include "G4OpProcessSubType.hh"
48#include "G4Poisson.hh"
52
53//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
54G4OpWLS::G4OpWLS(const G4String& processName, G4ProcessType type)
55 : G4VDiscreteProcess(processName, type)
56{
58 Initialise();
60 theIntegralTable = nullptr;
61
62 if(verboseLevel > 0)
63 G4cout << GetProcessName() << " is created " << G4endl;
64}
65
66//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
68{
70 {
72 delete theIntegralTable;
73 }
75}
76
77//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
79
80//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
82{
86}
87
88//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
90 const G4Step& aStep)
91{
92 std::vector<G4Track*> proposedSecondaries;
95
96 if(verboseLevel > 1)
97 {
98 G4cout << "\n** G4OpWLS: Photon absorbed! **" << G4endl;
99 }
100
101 G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
104 if(!MPT)
105 {
106 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
107 }
108 if(!MPT->GetProperty(kWLSCOMPONENT))
109 {
110 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
111 }
112
113 G4int NumPhotons = 1;
115 {
116 G4double MeanNumberOfPhotons = MPT->GetConstProperty(kWLSMEANNUMBERPHOTONS);
117 NumPhotons = G4int(G4Poisson(MeanNumberOfPhotons));
118 if(NumPhotons <= 0)
119 {
120 // return unchanged particle and no secondaries
122 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
123 }
124 }
125
126 // Retrieve the WLS Integral for this material
127 // new G4PhysicsFreeVector allocated to hold CII's
128 G4double primaryEnergy = aTrack.GetDynamicParticle()->GetKineticEnergy();
129 G4double WLSTime = 0.;
130 G4PhysicsFreeVector* WLSIntegral = nullptr;
131
132 WLSTime = MPT->GetConstProperty(kWLSTIMECONSTANT);
133 WLSIntegral = (G4PhysicsFreeVector*) ((*theIntegralTable)(
134 aTrack.GetMaterial()->GetIndex()));
135
136 // Max WLS Integral
137 G4double CIImax = WLSIntegral->GetMaxValue();
138 G4int NumberOfPhotons = NumPhotons;
139
140 for(G4int i = 0; i < NumPhotons; ++i)
141 {
142 G4double sampledEnergy;
143 // Make sure the energy of the secondary is less than that of the primary
144 for(G4int j = 1; j <= 100; ++j)
145 {
146 // Determine photon energy
147 G4double CIIvalue = G4UniformRand() * CIImax;
148 sampledEnergy = WLSIntegral->GetEnergy(CIIvalue);
149 if(sampledEnergy <= primaryEnergy)
150 break;
151 }
152 // If no such energy can be sampled, return one less secondary, or none
153 if(sampledEnergy > primaryEnergy)
154 {
155 if(verboseLevel > 1)
156 {
157 G4cout << " *** G4OpWLS: One less WLS photon will be returned ***"
158 << G4endl;
159 }
160 NumberOfPhotons--;
161 if(NumberOfPhotons == 0)
162 {
163 if(verboseLevel > 1)
164 {
165 G4cout
166 << " *** G4OpWLS: No WLS photon can be sampled for this primary ***"
167 << G4endl;
168 }
169 // return unchanged particle and no secondaries
171 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
172 }
173 continue;
174 }
175 else if(verboseLevel > 1)
176 {
177 G4cout << "G4OpWLS: Created photon with energy: " << sampledEnergy
178 << G4endl;
179 }
180
181 // Generate random photon direction
182 G4double cost = 1. - 2. * G4UniformRand();
183 G4double sint = std::sqrt((1. - cost) * (1. + cost));
184 G4double phi = twopi * G4UniformRand();
185 G4double sinp = std::sin(phi);
186 G4double cosp = std::cos(phi);
187 G4ParticleMomentum photonMomentum(sint * cosp, sint * sinp, cost);
188
189 G4ThreeVector photonPolarization(cost * cosp, cost * sinp, -sint);
190 G4ThreeVector perp = photonMomentum.cross(photonPolarization);
191
192 phi = twopi * G4UniformRand();
193 sinp = std::sin(phi);
194 cosp = std::cos(phi);
195 photonPolarization = (cosp * photonPolarization + sinp * perp).unit();
196
197 // Generate a new photon:
198 auto sec_dp =
200 sec_dp->SetPolarization(photonPolarization);
201 sec_dp->SetKineticEnergy(sampledEnergy);
202
203 G4double secTime = pPostStepPoint->GetGlobalTime() +
205 G4ThreeVector secPos = pPostStepPoint->GetPosition();
206 G4Track* secTrack = new G4Track(sec_dp, secTime, secPos);
207
208 secTrack->SetTouchableHandle(aTrack.GetTouchableHandle());
209 secTrack->SetParentID(aTrack.GetTrackID());
210
211 proposedSecondaries.push_back(secTrack);
212 }
213
214 aParticleChange.SetNumberOfSecondaries((G4int)proposedSecondaries.size());
215 for(auto sec : proposedSecondaries)
216 {
218 }
219 if(verboseLevel > 1)
220 {
221 G4cout << "\n Exiting from G4OpWLS::DoIt -- NumberOfSecondaries = "
223 }
224
225 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
226}
227
228//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
230{
232 {
234 delete theIntegralTable;
235 theIntegralTable = nullptr;
236 }
237
238 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
239 std::size_t numOfMaterials = G4Material::GetNumberOfMaterials();
240 theIntegralTable = new G4PhysicsTable(numOfMaterials);
241
242 // loop for materials
243 for(std::size_t i = 0; i < numOfMaterials; ++i)
244 {
245 auto physVector = new G4PhysicsFreeVector();
246
247 // Retrieve vector of WLS wavelength intensity for
248 // the material from the material's optical properties table.
250 (*materialTable)[i]->GetMaterialPropertiesTable();
251 if(MPT)
252 {
254 if(wlsVector)
255 {
256 // Retrieve the first intensity point in vector
257 // of (photon energy, intensity) pairs
258 G4double currentIN = (*wlsVector)[0];
259 if(currentIN >= 0.0)
260 {
261 // Create first (photon energy)
262 G4double currentPM = wlsVector->Energy(0);
263 G4double currentCII = 0.0;
264 physVector->InsertValues(currentPM, currentCII);
265
266 // Set previous values to current ones prior to loop
267 G4double prevPM = currentPM;
268 G4double prevCII = currentCII;
269 G4double prevIN = currentIN;
270
271 // loop over all (photon energy, intensity)
272 // pairs stored for this material
273 for(std::size_t j = 1; j < wlsVector->GetVectorLength(); ++j)
274 {
275 currentPM = wlsVector->Energy(j);
276 currentIN = (*wlsVector)[j];
277 currentCII =
278 prevCII + 0.5 * (currentPM - prevPM) * (prevIN + currentIN);
279
280 physVector->InsertValues(currentPM, currentCII);
281
282 prevPM = currentPM;
283 prevCII = currentCII;
284 prevIN = currentIN;
285 }
286 }
287 }
288 }
289 theIntegralTable->insertAt(i, physVector);
290 }
291}
292
293//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
296{
297 G4double thePhotonEnergy = aTrack.GetDynamicParticle()->GetTotalEnergy();
298 G4double attLength = DBL_MAX;
301
302 if(MPT)
303 {
305 if(attVector)
306 {
307 attLength = attVector->Value(thePhotonEnergy, idx_wls);
308 }
309 }
310 return attLength;
311}
312
313//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
315{
317 {
319 WLSTimeGeneratorProfile = nullptr;
320 }
321 if(name == "delta")
322 {
324 }
325 else if(name == "exponential")
326 {
328 new G4WLSTimeGeneratorProfileExponential("exponential");
329 }
330 else
331 {
332 G4Exception("G4OpWLS::UseTimeProfile", "em0202", FatalException,
333 "generator does not exist");
334 }
336}
337
338//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
340{
341 verboseLevel = verbose;
343}
@ FatalException
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:59
G4ForceCondition
@ kWLSMEANNUMBERPHOTONS
std::vector< G4Material * > G4MaterialTable
@ fOpWLS
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:50
G4ProcessType
@ fStopAndKill
double G4double
Definition: G4Types.hh:83
int G4int
Definition: G4Types.hh:85
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
#define G4UniformRand()
Definition: Randomize.hh:52
Hep3Vector cross(const Hep3Vector &) const
G4double GetKineticEnergy() const
G4double GetTotalEnergy() const
G4bool ConstPropertyExists(const G4String &key) const
G4double GetConstProperty(const G4String &key) const
G4MaterialPropertyVector * GetProperty(const char *key) const
G4MaterialPropertiesTable * GetMaterialPropertiesTable() const
Definition: G4Material.hh:251
static size_t GetNumberOfMaterials()
Definition: G4Material.cc:684
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:677
size_t GetIndex() const
Definition: G4Material.hh:255
virtual void BuildPhysicsTable(const G4ParticleDefinition &aParticleType) override
Definition: G4OpWLS.cc:229
virtual void Initialise()
Definition: G4OpWLS.cc:81
void SetVerboseLevel(G4int)
Definition: G4OpWLS.cc:339
virtual G4double GetMeanFreePath(const G4Track &aTrack, G4double, G4ForceCondition *) override
Definition: G4OpWLS.cc:294
G4PhysicsTable * theIntegralTable
Definition: G4OpWLS.hh:91
virtual G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep) override
Definition: G4OpWLS.cc:89
virtual ~G4OpWLS()
Definition: G4OpWLS.cc:67
G4VWLSTimeGeneratorProfile * WLSTimeGeneratorProfile
Definition: G4OpWLS.hh:90
G4OpWLS(const G4String &processName="OpWLS", G4ProcessType type=fOptical)
Definition: G4OpWLS.cc:54
virtual void UseTimeProfile(const G4String name)
Definition: G4OpWLS.cc:314
virtual void PreparePhysicsTable(const G4ParticleDefinition &) override
Definition: G4OpWLS.cc:78
static G4OpticalParameters * Instance()
G4String GetWLSTimeProfile() const
void SetWLSTimeProfile(const G4String &)
G4int GetWLSVerboseLevel() const
static G4OpticalPhoton * OpticalPhoton()
void AddSecondary(G4Track *aSecondary)
void Initialize(const G4Track &) override
void clearAndDestroy()
void insertAt(std::size_t, G4PhysicsVector *)
G4double GetEnergy(const G4double value) const
G4double GetMaxValue() const
G4double Energy(const std::size_t index) const
G4double Value(const G4double energy, std::size_t &lastidx) const
std::size_t GetVectorLength() const
G4double GetGlobalTime() const
const G4ThreeVector & GetPosition() const
Definition: G4Step.hh:62
G4StepPoint * GetPostStepPoint() const
G4int GetTrackID() const
void SetTouchableHandle(const G4TouchableHandle &apValue)
G4Material * GetMaterial() const
const G4DynamicParticle * GetDynamicParticle() const
const G4TouchableHandle & GetTouchableHandle() const
void SetParentID(const G4int aValue)
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
void ProposeTrackStatus(G4TrackStatus status)
G4int GetNumberOfSecondaries() const
void SetNumberOfSecondaries(G4int totSecondaries)
G4ParticleChange aParticleChange
Definition: G4VProcess.hh:331
G4int verboseLevel
Definition: G4VProcess.hh:360
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:410
const G4String & GetProcessName() const
Definition: G4VProcess.hh:386
virtual G4double GenerateTime(const G4double time_constant)=0
#define DBL_MAX
Definition: templates.hh:62