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

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

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