| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 //
28 ////////////////////////////////////////////////////////////////////////
29 // Optical Photon WaveLength Shifting (WLS) Class Implementation
30 ////////////////////////////////////////////////////////////////////////
31 //
32 // File: G4OpWLS2.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 "G4OpWLS2.hh"
44 #include "G4ios.hh"
45 #include "G4PhysicalConstants.hh"
46 #include "G4SystemOfUnits.hh"
47 #include "G4OpProcessSubType.hh"
48 #include "G4Poisson.hh"
49 #include "G4OpticalParameters.hh"
50 #include "G4WLSTimeGeneratorProfileDelta.hh"
51 #include "G4WLSTimeGeneratorProfileExponential.hh"
52
53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
54 G4OpWLS2::G4OpWLS2(const G4String& processName, G4ProcessType type)
55 : G4VDiscreteProcess(processName, type)
56 {
57 WLSTimeGeneratorProfile = nullptr;
58 Initialise();
59 SetProcessSubType(fOpWLS2);
60 theIntegralTable = nullptr;
61
62 if(verboseLevel > 0)
63 G4cout << GetProcessName() << " is created " << G4endl;
64 }
65
66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
67 G4OpWLS2::~G4OpWLS2()
68 {
69 if(theIntegralTable)
70 {
71 theIntegralTable->clearAndDestroy();
72 delete theIntegralTable;
73 }
74 delete WLSTimeGeneratorProfile;
75 }
76
77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
78 void G4OpWLS2::PreparePhysicsTable(const G4ParticleDefinition&)
79 {
80 Initialise();
81 }
82
83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
84 void G4OpWLS2::Initialise()
85 {
86 G4OpticalParameters* params = G4OpticalParameters::Instance();
87 SetVerboseLevel(params->GetWLS2VerboseLevel());
88 UseTimeProfile(params->GetWLS2TimeProfile());
89 }
90
91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
92 G4VParticleChange* G4OpWLS2::PostStepDoIt(const G4Track& aTrack,
93 const G4Step& aStep)
94 {
95 std::vector<G4Track*> proposedSecondaries;
96 aParticleChange.Initialize(aTrack);
97 aParticleChange.ProposeTrackStatus(fStopAndKill);
98
99 if(verboseLevel > 1)
100 {
101 G4cout << "\n** G4OpWLS2: Photon absorbed! **" << G4endl;
102 }
103
104 G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
105 G4MaterialPropertiesTable* MPT =
106 aTrack.GetMaterial()->GetMaterialPropertiesTable();
107 if(!MPT)
108 {
109 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
110 }
111 if(!MPT->GetProperty(kWLSCOMPONENT2))
112 {
113 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
114 }
115
116 G4int NumPhotons = 1;
117 if(MPT->ConstPropertyExists(kWLSMEANNUMBERPHOTONS2))
118 {
119 G4double MeanNumberOfPhotons =
120 MPT->GetConstProperty(kWLSMEANNUMBERPHOTONS2);
121 NumPhotons = G4int(G4Poisson(MeanNumberOfPhotons));
122 if(NumPhotons <= 0)
123 {
124 // return unchanged particle and no secondaries
125 aParticleChange.SetNumberOfSecondaries(0);
126 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
127 }
128 }
129
130 // Retrieve the WLS Integral for this material
131 // new G4PhysicsFreeVector allocated to hold CII's
132 G4double primaryEnergy = aTrack.GetDynamicParticle()->GetKineticEnergy();
133 G4double WLSTime = 0.;
134 G4PhysicsFreeVector* WLSIntegral = nullptr;
135
136 WLSTime = MPT->GetConstProperty(kWLSTIMECONSTANT2);
137 WLSIntegral = (G4PhysicsFreeVector*) ((*theIntegralTable)(
138 aTrack.GetMaterial()->GetIndex()));
139
140 // Max WLS Integral
141 G4double CIImax = WLSIntegral->GetMaxValue();
142 G4int NumberOfPhotons = NumPhotons;
143
144 for(G4int i = 0; i < NumPhotons; ++i)
145 {
146 G4double sampledEnergy;
147 // Make sure the energy of the secondary is less than that of the primary
148 for(G4int j = 1; j <= 100; ++j)
149 {
150 // Determine photon energy
151 G4double CIIvalue = G4UniformRand() * CIImax;
152 sampledEnergy = WLSIntegral->GetEnergy(CIIvalue);
153 if(sampledEnergy <= primaryEnergy)
154 break;
155 }
156 // If no such energy can be sampled, return one less secondary, or none
157 if(sampledEnergy > primaryEnergy)
158 {
159 if(verboseLevel > 1)
160 {
161 G4cout << " *** G4OpWLS2: One less WLS2 photon will be returned ***"
162 << G4endl;
163 }
164 NumberOfPhotons--;
165 if(NumberOfPhotons == 0)
166 {
167 if(verboseLevel > 1)
168 {
169 G4cout << " *** G4OpWLS2: No WLS2 photon can be sampled for this "
170 "primary ***"
171 << G4endl;
172 }
173 // return unchanged particle and no secondaries
174 aParticleChange.SetNumberOfSecondaries(0);
175 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
176 }
177 continue;
178 }
179 else if(verboseLevel > 1)
180 {
181 G4cout << "G4OpWLS2: Created photon with energy: " << sampledEnergy
182 << G4endl;
183 }
184
185 // Generate random photon direction
186 G4double cost = 1. - 2. * G4UniformRand();
187 G4double sint = std::sqrt((1. - cost) * (1. + cost));
188 G4double phi = twopi * G4UniformRand();
189 G4double sinp = std::sin(phi);
190 G4double cosp = std::cos(phi);
191 G4ParticleMomentum photonMomentum(sint * cosp, sint * sinp, cost);
192
193 G4ThreeVector photonPolarization(cost * cosp, cost * sinp, -sint);
194 G4ThreeVector perp = photonMomentum.cross(photonPolarization);
195
196 phi = twopi * G4UniformRand();
197 sinp = std::sin(phi);
198 cosp = std::cos(phi);
199 photonPolarization = (cosp * photonPolarization + sinp * perp).unit();
200
201 // Generate a new photon:
202 auto sec_dp =
203 new G4DynamicParticle(G4OpticalPhoton::OpticalPhoton(), photonMomentum);
204 sec_dp->SetPolarization(photonPolarization);
205 sec_dp->SetKineticEnergy(sampledEnergy);
206
207 G4double secTime = pPostStepPoint->GetGlobalTime() +
208 WLSTimeGeneratorProfile->GenerateTime(WLSTime);
209 G4ThreeVector secPos = pPostStepPoint->GetPosition();
210 G4Track* secTrack = new G4Track(sec_dp, secTime, secPos);
211
212 secTrack->SetTouchableHandle(aTrack.GetTouchableHandle());
213 secTrack->SetParentID(aTrack.GetTrackID());
214
215 proposedSecondaries.push_back(secTrack);
216 }
217
218 aParticleChange.SetNumberOfSecondaries((G4int)proposedSecondaries.size());
219 for(auto sec : proposedSecondaries)
220 {
221 aParticleChange.AddSecondary(sec);
222 }
223 if(verboseLevel > 1)
224 {
225 G4cout << "\n Exiting from G4OpWLS2::DoIt -- NumberOfSecondaries = "
226 << aParticleChange.GetNumberOfSecondaries() << G4endl;
227 }
228
229 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
230 }
231
232 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
233 void G4OpWLS2::BuildPhysicsTable(const G4ParticleDefinition&)
234 {
235 if(theIntegralTable)
236 {
237 theIntegralTable->clearAndDestroy();
238 delete theIntegralTable;
239 theIntegralTable = nullptr;
240 }
241
242 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
243 std::size_t numOfMaterials = G4Material::GetNumberOfMaterials();
244 theIntegralTable = new G4PhysicsTable(numOfMaterials);
245
246 // loop for materials
247 for(std::size_t i = 0; i < numOfMaterials; ++i)
248 {
249 auto physVector = new G4PhysicsFreeVector();
250
251 // Retrieve vector of WLS2 wavelength intensity for
252 // the material from the material's optical properties table.
253 G4MaterialPropertiesTable* MPT =
254 (*materialTable)[i]->GetMaterialPropertiesTable();
255 if(MPT)
256 {
257 G4MaterialPropertyVector* wlsVector = MPT->GetProperty(kWLSCOMPONENT2);
258 if(wlsVector)
259 {
260 // Retrieve the first intensity point in vector
261 // of (photon energy, intensity) pairs
262 G4double currentIN = (*wlsVector)[0];
263 if(currentIN >= 0.0)
264 {
265 // Create first (photon energy)
266 G4double currentPM = wlsVector->Energy(0);
267 G4double currentCII = 0.0;
268 physVector->InsertValues(currentPM, currentCII);
269
270 // Set previous values to current ones prior to loop
271 G4double prevPM = currentPM;
272 G4double prevCII = currentCII;
273 G4double prevIN = currentIN;
274
275 // loop over all (photon energy, intensity)
276 // pairs stored for this material
277 for(std::size_t j = 1; j < wlsVector->GetVectorLength(); ++j)
278 {
279 currentPM = wlsVector->Energy(j);
280 currentIN = (*wlsVector)[j];
281 currentCII =
282 prevCII + 0.5 * (currentPM - prevPM) * (prevIN + currentIN);
283
284 physVector->InsertValues(currentPM, currentCII);
285
286 prevPM = currentPM;
287 prevCII = currentCII;
288 prevIN = currentIN;
289 }
290 }
291 }
292 }
293 theIntegralTable->insertAt(i, physVector);
294 }
295 }
296
297 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
298 G4double G4OpWLS2::GetMeanFreePath(const G4Track& aTrack, G4double,
299 G4ForceCondition*)
300 {
301 G4double thePhotonEnergy = aTrack.GetDynamicParticle()->GetTotalEnergy();
302 G4double attLength = DBL_MAX;
303 G4MaterialPropertiesTable* MPT =
304 aTrack.GetMaterial()->GetMaterialPropertiesTable();
305
306 if(MPT)
307 {
308 G4MaterialPropertyVector* attVector = MPT->GetProperty(kWLSABSLENGTH2);
309 if(attVector)
310 {
311 attLength = attVector->Value(thePhotonEnergy, idx_wls2);
312 }
313 }
314 return attLength;
315 }
316
317 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
318 void G4OpWLS2::UseTimeProfile(const G4String name)
319 {
320 if(WLSTimeGeneratorProfile)
321 {
322 delete WLSTimeGeneratorProfile;
323 WLSTimeGeneratorProfile = nullptr;
324 }
325 if(name == "delta")
326 {
327 WLSTimeGeneratorProfile = new G4WLSTimeGeneratorProfileDelta("delta");
328 }
329 else if(name == "exponential")
330 {
331 WLSTimeGeneratorProfile =
332 new G4WLSTimeGeneratorProfileExponential("exponential");
333 }
334 else
335 {
336 G4Exception("G4OpWLS::UseTimeProfile", "em0202", FatalException,
337 "generator does not exist");
338 }
339 G4OpticalParameters::Instance()->SetWLS2TimeProfile(name);
340 }
341
342 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
343 void G4OpWLS2::SetVerboseLevel(G4int verbose)
344 {
345 verboseLevel = verbose;
346 G4OpticalParameters::Instance()->SetWLS2VerboseLevel(verboseLevel);
347 }
348