| 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: G4OpWLS2.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 "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
52
53 //....oooOO0OOooo........oooOO0OOooo........oo
54 G4OpWLS2::G4OpWLS2(const G4String& processName
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
64 }
65
66 //....oooOO0OOooo........oooOO0OOooo........oo
67 G4OpWLS2::~G4OpWLS2()
68 {
69 if(theIntegralTable)
70 {
71 theIntegralTable->clearAndDestroy();
72 delete theIntegralTable;
73 }
74 delete WLSTimeGeneratorProfile;
75 }
76
77 //....oooOO0OOooo........oooOO0OOooo........oo
78 void G4OpWLS2::PreparePhysicsTable(const G4Par
79 {
80 Initialise();
81 }
82
83 //....oooOO0OOooo........oooOO0OOooo........oo
84 void G4OpWLS2::Initialise()
85 {
86 G4OpticalParameters* params = G4OpticalParam
87 SetVerboseLevel(params->GetWLS2VerboseLevel(
88 UseTimeProfile(params->GetWLS2TimeProfile())
89 }
90
91 //....oooOO0OOooo........oooOO0OOooo........oo
92 G4VParticleChange* G4OpWLS2::PostStepDoIt(cons
93 cons
94 {
95 std::vector<G4Track*> proposedSecondaries;
96 aParticleChange.Initialize(aTrack);
97 aParticleChange.ProposeTrackStatus(fStopAndK
98
99 if(verboseLevel > 1)
100 {
101 G4cout << "\n** G4OpWLS2: Photon absorbed!
102 }
103
104 G4StepPoint* pPostStepPoint = aStep.GetPostS
105 G4MaterialPropertiesTable* MPT =
106 aTrack.GetMaterial()->GetMaterialPropertie
107 if(!MPT)
108 {
109 return G4VDiscreteProcess::PostStepDoIt(aT
110 }
111 if(!MPT->GetProperty(kWLSCOMPONENT2))
112 {
113 return G4VDiscreteProcess::PostStepDoIt(aT
114 }
115
116 G4int NumPhotons = 1;
117 if(MPT->ConstPropertyExists(kWLSMEANNUMBERPH
118 {
119 G4double MeanNumberOfPhotons =
120 MPT->GetConstProperty(kWLSMEANNUMBERPHOT
121 NumPhotons = G4int(G4Poisson(MeanNumberOfP
122 if(NumPhotons <= 0)
123 {
124 // return unchanged particle and no seco
125 aParticleChange.SetNumberOfSecondaries(0
126 return G4VDiscreteProcess::PostStepDoIt(
127 }
128 }
129
130 // Retrieve the WLS Integral for this materi
131 // new G4PhysicsFreeVector allocated to hold
132 G4double primaryEnergy = aTrack.GetDynamicPa
133 G4double WLSTime = 0.;
134 G4PhysicsFreeVector* WLSIntegral = nullptr;
135
136 WLSTime = MPT->GetConstProperty(kWLSTIME
137 WLSIntegral = (G4PhysicsFreeVector*) ((*theI
138 aTrack.GetMaterial()->GetIndex()));
139
140 // Max WLS Integral
141 G4double CIImax = WLSIntegral->GetMaxV
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 i
148 for(G4int j = 1; j <= 100; ++j)
149 {
150 // Determine photon energy
151 G4double CIIvalue = G4UniformRand() * CI
152 sampledEnergy = WLSIntegral->GetEner
153 if(sampledEnergy <= primaryEnergy)
154 break;
155 }
156 // If no such energy can be sampled, retur
157 if(sampledEnergy > primaryEnergy)
158 {
159 if(verboseLevel > 1)
160 {
161 G4cout << " *** G4OpWLS2: One less WLS
162 << G4endl;
163 }
164 NumberOfPhotons--;
165 if(NumberOfPhotons == 0)
166 {
167 if(verboseLevel > 1)
168 {
169 G4cout << " *** G4OpWLS2: No WLS2 ph
170 "primary ***"
171 << G4endl;
172 }
173 // return unchanged particle and no se
174 aParticleChange.SetNumberOfSecondaries
175 return G4VDiscreteProcess::PostStepDoI
176 }
177 continue;
178 }
179 else if(verboseLevel > 1)
180 {
181 G4cout << "G4OpWLS2: Created photon with
182 << G4endl;
183 }
184
185 // Generate random photon direction
186 G4double cost = 1. - 2. * G4UniformRand();
187 G4double sint = std::sqrt((1. - cost) * (1
188 G4double phi = twopi * G4UniformRand();
189 G4double sinp = std::sin(phi);
190 G4double cosp = std::cos(phi);
191 G4ParticleMomentum photonMomentum(sint * c
192
193 G4ThreeVector photonPolarization(cost * co
194 G4ThreeVector perp = photonMomentum.cross(
195
196 phi = twopi * G4UniformRand
197 sinp = std::sin(phi);
198 cosp = std::cos(phi);
199 photonPolarization = (cosp * photonPolariz
200
201 // Generate a new photon:
202 auto sec_dp =
203 new G4DynamicParticle(G4OpticalPhoton::O
204 sec_dp->SetPolarization(photonPolarization
205 sec_dp->SetKineticEnergy(sampledEnergy);
206
207 G4double secTime = pPostStepPoint->GetGlob
208 WLSTimeGeneratorProfile
209 G4ThreeVector secPos = pPostStepPoint->Get
210 G4Track* secTrack = new G4Track(sec_dp,
211
212 secTrack->SetTouchableHandle(aTrack.GetTou
213 secTrack->SetParentID(aTrack.GetTrackID())
214
215 proposedSecondaries.push_back(secTrack);
216 }
217
218 aParticleChange.SetNumberOfSecondaries((G4in
219 for(auto sec : proposedSecondaries)
220 {
221 aParticleChange.AddSecondary(sec);
222 }
223 if(verboseLevel > 1)
224 {
225 G4cout << "\n Exiting from G4OpWLS2::DoIt
226 << aParticleChange.GetNumberOfSecon
227 }
228
229 return G4VDiscreteProcess::PostStepDoIt(aTra
230 }
231
232 //....oooOO0OOooo........oooOO0OOooo........oo
233 void G4OpWLS2::BuildPhysicsTable(const G4Parti
234 {
235 if(theIntegralTable)
236 {
237 theIntegralTable->clearAndDestroy();
238 delete theIntegralTable;
239 theIntegralTable = nullptr;
240 }
241
242 const G4MaterialTable* materialTable = G4Mat
243 std::size_t numOfMaterials = G4Mat
244 theIntegralTable = new G
245
246 // loop for materials
247 for(std::size_t i = 0; i < numOfMaterials; +
248 {
249 auto physVector = new G4PhysicsFreeVector(
250
251 // Retrieve vector of WLS2 wavelength inte
252 // the material from the material's optica
253 G4MaterialPropertiesTable* MPT =
254 (*materialTable)[i]->GetMaterialProperti
255 if(MPT)
256 {
257 G4MaterialPropertyVector* wlsVector = MP
258 if(wlsVector)
259 {
260 // Retrieve the first intensity point
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->Ene
267 G4double currentCII = 0.0;
268 physVector->InsertValues(currentPM,
269
270 // Set previous values to current on
271 G4double prevPM = currentPM;
272 G4double prevCII = currentCII;
273 G4double prevIN = currentIN;
274
275 // loop over all (photon energy, int
276 // pairs stored for this material
277 for(std::size_t j = 1; j < wlsVector
278 {
279 currentPM = wlsVector->Energy(j);
280 currentIN = (*wlsVector)[j];
281 currentCII =
282 prevCII + 0.5 * (currentPM - pre
283
284 physVector->InsertValues(currentPM
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........oo
298 G4double G4OpWLS2::GetMeanFreePath(const G4Tra
299 G4ForceCond
300 {
301 G4double thePhotonEnergy = aTrack.GetDynamic
302 G4double attLength = DBL_MAX;
303 G4MaterialPropertiesTable* MPT =
304 aTrack.GetMaterial()->GetMaterialPropertie
305
306 if(MPT)
307 {
308 G4MaterialPropertyVector* attVector = MPT-
309 if(attVector)
310 {
311 attLength = attVector->Value(thePhotonEn
312 }
313 }
314 return attLength;
315 }
316
317 //....oooOO0OOooo........oooOO0OOooo........oo
318 void G4OpWLS2::UseTimeProfile(const G4String n
319 {
320 if(WLSTimeGeneratorProfile)
321 {
322 delete WLSTimeGeneratorProfile;
323 WLSTimeGeneratorProfile = nullptr;
324 }
325 if(name == "delta")
326 {
327 WLSTimeGeneratorProfile = new G4WLSTimeGen
328 }
329 else if(name == "exponential")
330 {
331 WLSTimeGeneratorProfile =
332 new G4WLSTimeGeneratorProfileExponential
333 }
334 else
335 {
336 G4Exception("G4OpWLS::UseTimeProfile", "em
337 "generator does not exist");
338 }
339 G4OpticalParameters::Instance()->SetWLS2Time
340 }
341
342 //....oooOO0OOooo........oooOO0OOooo........oo
343 void G4OpWLS2::SetVerboseLevel(G4int verbose)
344 {
345 verboseLevel = verbose;
346 G4OpticalParameters::Instance()->SetWLS2Verb
347 }
348