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

Last change on this file since 1244 was 1228, checked in by garnier, 14 years ago

update geant4.9.3 tag
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-03 $
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.