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

Last change on this file since 846 was 819, checked in by garnier, 16 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.