source: trunk/source/processes/optical/src/G4OpWLS.cc@ 1038

Last change on this file since 1038 was 1007, checked in by garnier, 17 years ago

update to geant4.9.2
File size: 11.6 KB
Line 
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// $Id: G4OpWLS.cc,v 1.13 2008/10/24 19:50:50 gum Exp $
28// GEANT4 tag $Name: geant4-09-02 $
29//
30////////////////////////////////////////////////////////////////////////
31// Optical Photon WaveLength Shifting (WLS) Class Implementation
32////////////////////////////////////////////////////////////////////////
33//
34// File: G4OpWLS.cc
35// Description: Discrete Process -- Wavelength Shifting of Optical Photons
36// Version: 1.0
37// Created: 2003-05-13
38// Author: John Paul Archambault
39// (Adaptation of G4Scintillation and G4OpAbsorption)
40// Updated: 2005-07-28 - add G4ProcessType to constructor
41// 2006-05-07 - add G4VWLSTimeGeneratorProfile
42// mail: gum@triumf.ca
43// jparcham@phys.ualberta.ca
44//
45////////////////////////////////////////////////////////////////////////
46
47#include "G4ios.hh"
48#include "G4OpProcessSubType.hh"
49
50#include "G4OpWLS.hh"
51#include "G4WLSTimeGeneratorProfileDelta.hh"
52#include "G4WLSTimeGeneratorProfileExponential.hh"
53
54/////////////////////////
55// Class Implementation
56/////////////////////////
57
58/////////////////
59// Constructors
60/////////////////
61
62G4OpWLS::G4OpWLS(const G4String& processName, G4ProcessType type)
63 : G4VDiscreteProcess(processName, type)
64{
65 SetProcessSubType(fOpWLS);
66
67 theIntegralTable = 0;
68
69 if (verboseLevel>0) {
70 G4cout << GetProcessName() << " is created " << G4endl;
71 }
72
73 WLSTimeGeneratorProfile =
74 new G4WLSTimeGeneratorProfileDelta("WLSTimeGeneratorProfileDelta");
75
76 BuildThePhysicsTable();
77}
78
79////////////////
80// Destructors
81////////////////
82
83G4OpWLS::~G4OpWLS()
84{
85 if (theIntegralTable != 0) {
86 theIntegralTable->clearAndDestroy();
87 delete theIntegralTable;
88 }
89 delete WLSTimeGeneratorProfile;
90}
91
92////////////
93// Methods
94////////////
95
96// PostStepDoIt
97// -------------
98//
99G4VParticleChange*
100G4OpWLS::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
101{
102 aParticleChange.Initialize(aTrack);
103
104 aParticleChange.ProposeTrackStatus(fStopAndKill);
105
106 if (verboseLevel>0) {
107 G4cout << "\n** Photon absorbed! **" << G4endl;
108 }
109
110 const G4Material* aMaterial = aTrack.GetMaterial();
111
112 G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
113
114 G4MaterialPropertiesTable* aMaterialPropertiesTable =
115 aMaterial->GetMaterialPropertiesTable();
116 if (!aMaterialPropertiesTable)
117 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
118
119 const G4MaterialPropertyVector* WLS_Intensity =
120 aMaterialPropertiesTable->GetProperty("WLSCOMPONENT");
121
122 if (!WLS_Intensity)
123 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
124
125 G4int NumPhotons = 1;
126
127 if (aMaterialPropertiesTable->ConstPropertyExists("WLSMEANNUMBERPHOTONS")) {
128
129 G4double MeanNumberOfPhotons = aMaterialPropertiesTable->
130 GetConstProperty("WLSMEANNUMBERPHOTONS");
131
132 NumPhotons = G4int(G4Poisson(MeanNumberOfPhotons));
133
134 if (NumPhotons <= 0) {
135
136 // return unchanged particle and no secondaries
137
138 aParticleChange.SetNumberOfSecondaries(0);
139
140 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
141
142 }
143
144 }
145
146 aParticleChange.SetNumberOfSecondaries(NumPhotons);
147
148 G4int materialIndex = aMaterial->GetIndex();
149
150 // Retrieve the WLS Integral for this material
151 // new G4PhysicsOrderedFreeVector allocated to hold CII's
152
153 G4double WLSTime = 0.*ns;
154 G4PhysicsOrderedFreeVector* WLSIntegral = 0;
155
156 WLSTime = aMaterialPropertiesTable->
157 GetConstProperty("WLSTIMECONSTANT");
158 WLSIntegral =
159 (G4PhysicsOrderedFreeVector*)((*theIntegralTable)(materialIndex));
160
161 // Max WLS Integral
162
163 G4double CIImax = WLSIntegral->GetMaxValue();
164
165 for (G4int i = 0; i < NumPhotons; i++) {
166
167 // Determine photon energy
168
169 G4double CIIvalue = G4UniformRand()*CIImax;
170 G4double sampledEnergy =
171 WLSIntegral->GetEnergy(CIIvalue);
172
173 if (verboseLevel>1) {
174 G4cout << "sampledEnergy = " << sampledEnergy << G4endl;
175 G4cout << "CIIvalue = " << CIIvalue << G4endl;
176 }
177
178 // Generate random photon direction
179
180 G4double cost = 1. - 2.*G4UniformRand();
181 G4double sint = std::sqrt((1.-cost)*(1.+cost));
182
183 G4double phi = twopi*G4UniformRand();
184 G4double sinp = std::sin(phi);
185 G4double cosp = std::cos(phi);
186
187 G4double px = sint*cosp;
188 G4double py = sint*sinp;
189 G4double pz = cost;
190
191 // Create photon momentum direction vector
192
193 G4ParticleMomentum photonMomentum(px, py, pz);
194
195 // Determine polarization of new photon
196
197 G4double sx = cost*cosp;
198 G4double sy = cost*sinp;
199 G4double sz = -sint;
200
201 G4ThreeVector photonPolarization(sx, sy, sz);
202
203 G4ThreeVector perp = photonMomentum.cross(photonPolarization);
204
205 phi = twopi*G4UniformRand();
206 sinp = std::sin(phi);
207 cosp = std::cos(phi);
208
209 photonPolarization = cosp * photonPolarization + sinp * perp;
210
211 photonPolarization = photonPolarization.unit();
212
213 // Generate a new photon:
214
215 G4DynamicParticle* aWLSPhoton =
216 new G4DynamicParticle(G4OpticalPhoton::OpticalPhoton(),
217 photonMomentum);
218 aWLSPhoton->SetPolarization
219 (photonPolarization.x(),
220 photonPolarization.y(),
221 photonPolarization.z());
222
223 aWLSPhoton->SetKineticEnergy(sampledEnergy);
224
225 // Generate new G4Track object:
226
227 // Must give position of WLS optical photon
228
229 G4double TimeDelay = WLSTimeGeneratorProfile->GenerateTime(WLSTime);
230 G4double aSecondaryTime = (pPostStepPoint->GetGlobalTime()) + TimeDelay;
231
232 G4ThreeVector aSecondaryPosition = pPostStepPoint->GetPosition();
233
234 G4Track* aSecondaryTrack =
235 new G4Track(aWLSPhoton,aSecondaryTime,aSecondaryPosition);
236
237 aSecondaryTrack->SetTouchableHandle(aTrack.GetTouchableHandle());
238 // aSecondaryTrack->SetTouchableHandle((G4VTouchable*)0);
239
240 aSecondaryTrack->SetParentID(aTrack.GetTrackID());
241
242 aParticleChange.AddSecondary(aSecondaryTrack);
243 }
244
245 if (verboseLevel>0) {
246 G4cout << "\n Exiting from G4OpWLS::DoIt -- NumberOfSecondaries = "
247 << aParticleChange.GetNumberOfSecondaries() << G4endl;
248 }
249
250 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
251}
252
253// BuildThePhysicsTable for the wavelength shifting process
254// --------------------------------------------------
255//
256
257void G4OpWLS::BuildThePhysicsTable()
258{
259 if (theIntegralTable) return;
260
261 const G4MaterialTable* theMaterialTable =
262 G4Material::GetMaterialTable();
263 G4int numOfMaterials = G4Material::GetNumberOfMaterials();
264
265 // create new physics table
266
267 if(!theIntegralTable)theIntegralTable = new G4PhysicsTable(numOfMaterials);
268
269 // loop for materials
270
271 for (G4int i=0 ; i < numOfMaterials; i++)
272 {
273 G4PhysicsOrderedFreeVector* aPhysicsOrderedFreeVector =
274 new G4PhysicsOrderedFreeVector();
275
276 // Retrieve vector of WLS wavelength intensity for
277 // the material from the material's optical properties table.
278
279 G4Material* aMaterial = (*theMaterialTable)[i];
280
281 G4MaterialPropertiesTable* aMaterialPropertiesTable =
282 aMaterial->GetMaterialPropertiesTable();
283
284 if (aMaterialPropertiesTable) {
285
286 G4MaterialPropertyVector* theWLSVector =
287 aMaterialPropertiesTable->GetProperty("WLSCOMPONENT");
288
289 if (theWLSVector) {
290
291 // Retrieve the first intensity point in vector
292 // of (photon energy, intensity) pairs
293
294 theWLSVector->ResetIterator();
295 ++(*theWLSVector); // advance to 1st entry
296
297 G4double currentIN = theWLSVector->
298 GetProperty();
299
300 if (currentIN >= 0.0) {
301
302 // Create first (photon energy)
303
304 G4double currentPM = theWLSVector->
305 GetPhotonEnergy();
306
307 G4double currentCII = 0.0;
308
309 aPhysicsOrderedFreeVector->
310 InsertValues(currentPM , currentCII);
311
312 // Set previous values to current ones prior to loop
313
314 G4double prevPM = currentPM;
315 G4double prevCII = currentCII;
316 G4double prevIN = currentIN;
317
318 // loop over all (photon energy, intensity)
319 // pairs stored for this material
320
321 while(++(*theWLSVector))
322 {
323 currentPM = theWLSVector->
324 GetPhotonEnergy();
325
326 currentIN=theWLSVector->
327 GetProperty();
328
329 currentCII = 0.5 * (prevIN + currentIN);
330
331 currentCII = prevCII +
332 (currentPM - prevPM) * currentCII;
333
334 aPhysicsOrderedFreeVector->
335 InsertValues(currentPM, currentCII);
336
337 prevPM = currentPM;
338 prevCII = currentCII;
339 prevIN = currentIN;
340 }
341 }
342 }
343 }
344 // The WLS integral for a given material
345 // will be inserted in the table according to the
346 // position of the material in the material table.
347
348 theIntegralTable->insertAt(i,aPhysicsOrderedFreeVector);
349 }
350}
351
352// GetMeanFreePath
353// ---------------
354//
355G4double G4OpWLS::GetMeanFreePath(const G4Track& aTrack,
356 G4double ,
357 G4ForceCondition* )
358{
359 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
360 const G4Material* aMaterial = aTrack.GetMaterial();
361
362 G4double thePhotonEnergy = aParticle->GetTotalEnergy();
363
364 G4MaterialPropertiesTable* aMaterialPropertyTable;
365 G4MaterialPropertyVector* AttenuationLengthVector;
366
367 G4double AttenuationLength = DBL_MAX;
368
369 aMaterialPropertyTable = aMaterial->GetMaterialPropertiesTable();
370
371 if ( aMaterialPropertyTable ) {
372 AttenuationLengthVector = aMaterialPropertyTable->
373 GetProperty("WLSABSLENGTH");
374 if ( AttenuationLengthVector ){
375 AttenuationLength = AttenuationLengthVector->
376 GetProperty (thePhotonEnergy);
377 }
378 else {
379 // G4cout << "No WLS absorption length specified" << G4endl;
380 }
381 }
382 else {
383 // G4cout << "No WLS absortion length specified" << G4endl;
384 }
385
386 return AttenuationLength;
387}
388
389void G4OpWLS::UseTimeProfile(const G4String name)
390{
391 if (name == "delta")
392 {
393 delete WLSTimeGeneratorProfile;
394 WLSTimeGeneratorProfile =
395 new G4WLSTimeGeneratorProfileDelta("delta");
396 }
397 else if (name == "exponential")
398 {
399 delete WLSTimeGeneratorProfile;
400 WLSTimeGeneratorProfile =
401 new G4WLSTimeGeneratorProfileExponential("exponential");
402 }
403 else
404 {
405 G4Exception("G4OpWLS::UseTimeProfile - generator does not exist");
406 }
407}
Note: See TracBrowser for help on using the repository browser.