| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 //
28 //////////////////////////////////////////////
29 // Optical Photon WaveLength Shifting (WLS) Cl
30 //////////////////////////////////////////////
31 //
32 // File: G4OpWLS.cc
33 // Description: Discrete Process -- Wavelength
34 // Version: 1.0
35 // Created: 2003-05-13
36 // Author: John Paul Archambault
37 // (Adaptation of G4Scintillation
38 // Updated: 2005-07-28 - add G4ProcessType
39 // 2006-05-07 - add G4VWLSTimeGen
40 //
41 //////////////////////////////////////////////
42
43 #include "G4OpWLS.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
52
53 //....oooOO0OOooo........oooOO0OOooo........oo
54 G4OpWLS::G4OpWLS(const G4String& processName,
55 : G4VDiscreteProcess(processName, type)
56 {
57 WLSTimeGeneratorProfile = nullptr;
58 Initialise();
59 SetProcessSubType(fOpWLS);
60 theIntegralTable = nullptr;
61
62 if(verboseLevel > 0)
63 G4cout << GetProcessName() << " is created
64 }
65
66 //....oooOO0OOooo........oooOO0OOooo........oo
67 G4OpWLS::~G4OpWLS()
68 {
69 if(theIntegralTable)
70 {
71 theIntegralTable->clearAndDestroy();
72 delete theIntegralTable;
73 }
74 delete WLSTimeGeneratorProfile;
75 }
76
77 //....oooOO0OOooo........oooOO0OOooo........oo
78 void G4OpWLS::PreparePhysicsTable(const G4Part
79
80 //....oooOO0OOooo........oooOO0OOooo........oo
81 void G4OpWLS::Initialise()
82 {
83 G4OpticalParameters* params = G4OpticalParam
84 SetVerboseLevel(params->GetWLSVerboseLevel()
85 UseTimeProfile(params->GetWLSTimeProfile());
86 }
87
88 //....oooOO0OOooo........oooOO0OOooo........oo
89 G4VParticleChange* G4OpWLS::PostStepDoIt(const
90 const
91 {
92 std::vector<G4Track*> proposedSecondaries;
93 aParticleChange.Initialize(aTrack);
94 aParticleChange.ProposeTrackStatus(fStopAndK
95
96 if(verboseLevel > 1)
97 {
98 G4cout << "\n** G4OpWLS: Photon absorbed!
99 }
100
101 G4StepPoint* pPostStepPoint = aStep.GetPostS
102 G4MaterialPropertiesTable* MPT =
103 aTrack.GetMaterial()->GetMaterialPropertie
104 if(!MPT)
105 {
106 return G4VDiscreteProcess::PostStepDoIt(aT
107 }
108 if(!MPT->GetProperty(kWLSCOMPONENT))
109 {
110 return G4VDiscreteProcess::PostStepDoIt(aT
111 }
112
113 G4int NumPhotons = 1;
114 if(MPT->ConstPropertyExists(kWLSMEANNUMBERPH
115 {
116 G4double MeanNumberOfPhotons = MPT->GetCon
117 NumPhotons = G4int(G4Poi
118 if(NumPhotons <= 0)
119 {
120 // return unchanged particle and no seco
121 aParticleChange.SetNumberOfSecondaries(0
122 return G4VDiscreteProcess::PostStepDoIt(
123 }
124 }
125
126 // Retrieve the WLS Integral for this materi
127 // new G4PhysicsFreeVector allocated to hold
128 G4double primaryEnergy = aTrack.GetDynamicPa
129 G4double WLSTime = 0.;
130 G4PhysicsFreeVector* WLSIntegral = nullptr;
131
132 WLSTime = MPT->GetConstProperty(kWLSTIME
133 WLSIntegral = (G4PhysicsFreeVector*) ((*theI
134 aTrack.GetMaterial()->GetIndex()));
135
136 // Max WLS Integral
137 G4double CIImax = WLSIntegral->GetMaxV
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 i
144 for(G4int j = 1; j <= 100; ++j)
145 {
146 // Determine photon energy
147 G4double CIIvalue = G4UniformRand() * CI
148 sampledEnergy = WLSIntegral->GetEner
149 if(sampledEnergy <= primaryEnergy)
150 break;
151 }
152 // If no such energy can be sampled, retur
153 if(sampledEnergy > primaryEnergy)
154 {
155 if(verboseLevel > 1)
156 {
157 G4cout << " *** G4OpWLS: One less WLS
158 << G4endl;
159 }
160 NumberOfPhotons--;
161 if(NumberOfPhotons == 0)
162 {
163 if(verboseLevel > 1)
164 {
165 G4cout
166 << " *** G4OpWLS: No WLS photon ca
167 << G4endl;
168 }
169 // return unchanged particle and no se
170 aParticleChange.SetNumberOfSecondaries
171 return G4VDiscreteProcess::PostStepDoI
172 }
173 continue;
174 }
175 else if(verboseLevel > 1)
176 {
177 G4cout << "G4OpWLS: Created photon with
178 << G4endl;
179 }
180
181 // Generate random photon direction
182 G4double cost = 1. - 2. * G4UniformRand();
183 G4double sint = std::sqrt((1. - cost) * (1
184 G4double phi = twopi * G4UniformRand();
185 G4double sinp = std::sin(phi);
186 G4double cosp = std::cos(phi);
187 G4ParticleMomentum photonMomentum(sint * c
188
189 G4ThreeVector photonPolarization(cost * co
190 G4ThreeVector perp = photonMomentum.cross(
191
192 phi = twopi * G4UniformRand
193 sinp = std::sin(phi);
194 cosp = std::cos(phi);
195 photonPolarization = (cosp * photonPolariz
196
197 // Generate a new photon:
198 auto sec_dp =
199 new G4DynamicParticle(G4OpticalPhoton::O
200 sec_dp->SetPolarization(photonPolarization
201 sec_dp->SetKineticEnergy(sampledEnergy);
202
203 G4double secTime = pPostStepPoint->GetGlob
204 WLSTimeGeneratorProfile
205 G4ThreeVector secPos = pPostStepPoint->Get
206 G4Track* secTrack = new G4Track(sec_dp,
207
208 secTrack->SetTouchableHandle(aTrack.GetTou
209 secTrack->SetParentID(aTrack.GetTrackID())
210
211 proposedSecondaries.push_back(secTrack);
212 }
213
214 aParticleChange.SetNumberOfSecondaries((G4in
215 for(auto sec : proposedSecondaries)
216 {
217 aParticleChange.AddSecondary(sec);
218 }
219 if(verboseLevel > 1)
220 {
221 G4cout << "\n Exiting from G4OpWLS::DoIt -
222 << aParticleChange.GetNumberOfSecon
223 }
224
225 return G4VDiscreteProcess::PostStepDoIt(aTra
226 }
227
228 //....oooOO0OOooo........oooOO0OOooo........oo
229 void G4OpWLS::BuildPhysicsTable(const G4Partic
230 {
231 if(theIntegralTable)
232 {
233 theIntegralTable->clearAndDestroy();
234 delete theIntegralTable;
235 theIntegralTable = nullptr;
236 }
237
238 const G4MaterialTable* materialTable = G4Mat
239 std::size_t numOfMaterials = G4Mat
240 theIntegralTable = new G
241
242 // loop for materials
243 for(std::size_t i = 0; i < numOfMaterials; +
244 {
245 auto physVector = new G4PhysicsFreeVector(
246
247 // Retrieve vector of WLS wavelength inten
248 // the material from the material's optica
249 G4MaterialPropertiesTable* MPT =
250 (*materialTable)[i]->GetMaterialProperti
251 if(MPT)
252 {
253 G4MaterialPropertyVector* wlsVector = MP
254 if(wlsVector)
255 {
256 // Retrieve the first intensity point
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->Ene
263 G4double currentCII = 0.0;
264 physVector->InsertValues(currentPM,
265
266 // Set previous values to current on
267 G4double prevPM = currentPM;
268 G4double prevCII = currentCII;
269 G4double prevIN = currentIN;
270
271 // loop over all (photon energy, int
272 // pairs stored for this material
273 for(std::size_t j = 1; j < wlsVector
274 {
275 currentPM = wlsVector->Energy(j);
276 currentIN = (*wlsVector)[j];
277 currentCII =
278 prevCII + 0.5 * (currentPM - pre
279
280 physVector->InsertValues(currentPM
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........oo
294 G4double G4OpWLS::GetMeanFreePath(const G4Trac
295 G4ForceCondi
296 {
297 G4double thePhotonEnergy = aTrack.GetDynamic
298 G4double attLength = DBL_MAX;
299 G4MaterialPropertiesTable* MPT =
300 aTrack.GetMaterial()->GetMaterialPropertie
301
302 if(MPT)
303 {
304 G4MaterialPropertyVector* attVector = MPT-
305 if(attVector)
306 {
307 attLength = attVector->Value(thePhotonEn
308 }
309 }
310 return attLength;
311 }
312
313 //....oooOO0OOooo........oooOO0OOooo........oo
314 void G4OpWLS::UseTimeProfile(const G4String na
315 {
316 if(WLSTimeGeneratorProfile)
317 {
318 delete WLSTimeGeneratorProfile;
319 WLSTimeGeneratorProfile = nullptr;
320 }
321 if(name == "delta")
322 {
323 WLSTimeGeneratorProfile = new G4WLSTimeGen
324 }
325 else if(name == "exponential")
326 {
327 WLSTimeGeneratorProfile =
328 new G4WLSTimeGeneratorProfileExponential
329 }
330 else
331 {
332 G4Exception("G4OpWLS::UseTimeProfile", "em
333 "generator does not exist");
334 }
335 G4OpticalParameters::Instance()->SetWLSTimeP
336 }
337
338 //....oooOO0OOooo........oooOO0OOooo........oo
339 void G4OpWLS::SetVerboseLevel(G4int verbose)
340 {
341 verboseLevel = verbose;
342 G4OpticalParameters::Instance()->SetWLSVerbo
343 }
344