Geant4  10.00.p01
G4OpWLS.cc
Go to the documentation of this file.
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 79221 2014-02-20 14:58:02Z gcosmo $
28 //
30 // Optical Photon WaveLength Shifting (WLS) Class Implementation
32 //
33 // File: G4OpWLS.cc
34 // Description: Discrete Process -- Wavelength Shifting of Optical Photons
35 // Version: 1.0
36 // Created: 2003-05-13
37 // Author: John Paul Archambault
38 // (Adaptation of G4Scintillation and G4OpAbsorption)
39 // Updated: 2005-07-28 - add G4ProcessType to constructor
40 // 2006-05-07 - add G4VWLSTimeGeneratorProfile
41 // mail: gum@triumf.ca
42 // jparcham@phys.ualberta.ca
43 //
45 
46 #include "G4OpWLS.hh"
47 
48 #include "G4ios.hh"
49 #include "G4PhysicalConstants.hh"
50 #include "G4SystemOfUnits.hh"
51 #include "G4OpProcessSubType.hh"
52 
55 
57 // Class Implementation
59 
61  // static data members
63 
65 
67 // Constructors
69 
70 G4OpWLS::G4OpWLS(const G4String& processName, G4ProcessType type)
71  : G4VDiscreteProcess(processName, type)
72 {
74 
75  theIntegralTable = NULL;
77  { WLSTimeGeneratorProfile = new G4WLSTimeGeneratorProfileDelta("WLSTimeGeneratorProfileDelta"); }
78 
79  if (verboseLevel>0) {
80  G4cout << GetProcessName() << " is created " << G4endl;
81  }
82 }
83 
85 // Destructors
87 
89 {
90  if (theIntegralTable != 0) {
92  delete theIntegralTable;
93  }
94 }
95 
97 // Methods
99 
101 {
103 }
104 
105 // PostStepDoIt
106 // -------------
107 //
109 G4OpWLS::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
110 {
111  aParticleChange.Initialize(aTrack);
112 
114 
115  if (verboseLevel>0) {
116  G4cout << "\n** Photon absorbed! **" << G4endl;
117  }
118 
119  const G4Material* aMaterial = aTrack.GetMaterial();
120 
121  G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
122 
123  G4MaterialPropertiesTable* aMaterialPropertiesTable =
124  aMaterial->GetMaterialPropertiesTable();
125  if (!aMaterialPropertiesTable)
126  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
127 
128  const G4MaterialPropertyVector* WLS_Intensity =
129  aMaterialPropertiesTable->GetProperty("WLSCOMPONENT");
130 
131  if (!WLS_Intensity)
132  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
133 
134  G4int NumPhotons = 1;
135 
136  if (aMaterialPropertiesTable->ConstPropertyExists("WLSMEANNUMBERPHOTONS")) {
137 
138  G4double MeanNumberOfPhotons = aMaterialPropertiesTable->
139  GetConstProperty("WLSMEANNUMBERPHOTONS");
140 
141  NumPhotons = G4int(G4Poisson(MeanNumberOfPhotons));
142 
143  if (NumPhotons <= 0) {
144 
145  // return unchanged particle and no secondaries
146 
148 
149  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
150 
151  }
152 
153  }
154 
156 
157  G4double primaryEnergy = aTrack.GetDynamicParticle()->GetKineticEnergy();
158 
159  G4int materialIndex = aMaterial->GetIndex();
160 
161  // Retrieve the WLS Integral for this material
162  // new G4PhysicsOrderedFreeVector allocated to hold CII's
163 
164  G4double WLSTime = 0.*ns;
165  G4PhysicsOrderedFreeVector* WLSIntegral = 0;
166 
167  WLSTime = aMaterialPropertiesTable->
168  GetConstProperty("WLSTIMECONSTANT");
169  WLSIntegral =
170  (G4PhysicsOrderedFreeVector*)((*theIntegralTable)(materialIndex));
171 
172  // Max WLS Integral
173 
174  G4double CIImax = WLSIntegral->GetMaxValue();
175 
176  G4int NumberOfPhotons = NumPhotons;
177 
178  for (G4int i = 0; i < NumPhotons; i++) {
179 
180  G4double sampledEnergy;
181 
182  // Make sure the energy of the secondary is less than that of the primary
183 
184  for (G4int j = 1; j <= 100; j++) {
185 
186  // Determine photon energy
187 
188  G4double CIIvalue = G4UniformRand()*CIImax;
189  sampledEnergy = WLSIntegral->GetEnergy(CIIvalue);
190 
191  if (verboseLevel>1) {
192  G4cout << "sampledEnergy = " << sampledEnergy << G4endl;
193  G4cout << "CIIvalue = " << CIIvalue << G4endl;
194  }
195 
196  if (sampledEnergy <= primaryEnergy) break;
197  }
198 
199  // If no such energy can be sampled, return one less secondary, or none
200 
201  if (sampledEnergy > primaryEnergy) {
202  if (verboseLevel>1)
203  G4cout << " *** One less WLS photon will be returned ***" << G4endl;
204  NumberOfPhotons--;
205  aParticleChange.SetNumberOfSecondaries(NumberOfPhotons);
206  if (NumberOfPhotons == 0) {
207  if (verboseLevel>1)
208  G4cout << " *** No WLS photon can be sampled for this primary ***"
209  << G4endl;
210  // return unchanged particle and no secondaries
211  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
212  }
213  continue;
214  }
215 
216  // Generate random photon direction
217 
218  G4double cost = 1. - 2.*G4UniformRand();
219  G4double sint = std::sqrt((1.-cost)*(1.+cost));
220 
221  G4double phi = twopi*G4UniformRand();
222  G4double sinp = std::sin(phi);
223  G4double cosp = std::cos(phi);
224 
225  G4double px = sint*cosp;
226  G4double py = sint*sinp;
227  G4double pz = cost;
228 
229  // Create photon momentum direction vector
230 
231  G4ParticleMomentum photonMomentum(px, py, pz);
232 
233  // Determine polarization of new photon
234 
235  G4double sx = cost*cosp;
236  G4double sy = cost*sinp;
237  G4double sz = -sint;
238 
239  G4ThreeVector photonPolarization(sx, sy, sz);
240 
241  G4ThreeVector perp = photonMomentum.cross(photonPolarization);
242 
243  phi = twopi*G4UniformRand();
244  sinp = std::sin(phi);
245  cosp = std::cos(phi);
246 
247  photonPolarization = cosp * photonPolarization + sinp * perp;
248 
249  photonPolarization = photonPolarization.unit();
250 
251  // Generate a new photon:
252 
253  G4DynamicParticle* aWLSPhoton =
255  photonMomentum);
256  aWLSPhoton->SetPolarization
257  (photonPolarization.x(),
258  photonPolarization.y(),
259  photonPolarization.z());
260 
261  aWLSPhoton->SetKineticEnergy(sampledEnergy);
262 
263  // Generate new G4Track object:
264 
265  // Must give position of WLS optical photon
266 
267  G4double TimeDelay = WLSTimeGeneratorProfile->GenerateTime(WLSTime);
268  G4double aSecondaryTime = (pPostStepPoint->GetGlobalTime()) + TimeDelay;
269 
270  G4ThreeVector aSecondaryPosition = pPostStepPoint->GetPosition();
271 
272  G4Track* aSecondaryTrack =
273  new G4Track(aWLSPhoton,aSecondaryTime,aSecondaryPosition);
274 
275  aSecondaryTrack->SetTouchableHandle(aTrack.GetTouchableHandle());
276  // aSecondaryTrack->SetTouchableHandle((G4VTouchable*)0);
277 
278  aSecondaryTrack->SetParentID(aTrack.GetTrackID());
279 
280  aParticleChange.AddSecondary(aSecondaryTrack);
281  }
282 
283  if (verboseLevel>0) {
284  G4cout << "\n Exiting from G4OpWLS::DoIt -- NumberOfSecondaries = "
286  }
287 
288  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
289 }
290 
291 // BuildThePhysicsTable for the wavelength shifting process
292 // --------------------------------------------------
293 //
294 
296 {
297  if (theIntegralTable) return;
298 
299  const G4MaterialTable* theMaterialTable =
301  G4int numOfMaterials = G4Material::GetNumberOfMaterials();
302 
303  // create new physics table
304 
305  if(!theIntegralTable)theIntegralTable = new G4PhysicsTable(numOfMaterials);
306 
307  // loop for materials
308 
309  for (G4int i=0 ; i < numOfMaterials; i++)
310  {
311  G4PhysicsOrderedFreeVector* aPhysicsOrderedFreeVector =
313 
314  // Retrieve vector of WLS wavelength intensity for
315  // the material from the material's optical properties table.
316 
317  G4Material* aMaterial = (*theMaterialTable)[i];
318 
319  G4MaterialPropertiesTable* aMaterialPropertiesTable =
320  aMaterial->GetMaterialPropertiesTable();
321 
322  if (aMaterialPropertiesTable) {
323 
324  G4MaterialPropertyVector* theWLSVector =
325  aMaterialPropertiesTable->GetProperty("WLSCOMPONENT");
326 
327  if (theWLSVector) {
328 
329  // Retrieve the first intensity point in vector
330  // of (photon energy, intensity) pairs
331 
332  G4double currentIN = (*theWLSVector)[0];
333 
334  if (currentIN >= 0.0) {
335 
336  // Create first (photon energy)
337 
338  G4double currentPM = theWLSVector->Energy(0);
339 
340  G4double currentCII = 0.0;
341 
342  aPhysicsOrderedFreeVector->
343  InsertValues(currentPM , currentCII);
344 
345  // Set previous values to current ones prior to loop
346 
347  G4double prevPM = currentPM;
348  G4double prevCII = currentCII;
349  G4double prevIN = currentIN;
350 
351  // loop over all (photon energy, intensity)
352  // pairs stored for this material
353 
354  for (size_t j = 1;
355  j < theWLSVector->GetVectorLength();
356  j++)
357  {
358  currentPM = theWLSVector->Energy(j);
359  currentIN = (*theWLSVector)[j];
360 
361  currentCII = 0.5 * (prevIN + currentIN);
362 
363  currentCII = prevCII +
364  (currentPM - prevPM) * currentCII;
365 
366  aPhysicsOrderedFreeVector->
367  InsertValues(currentPM, currentCII);
368 
369  prevPM = currentPM;
370  prevCII = currentCII;
371  prevIN = currentIN;
372  }
373  }
374  }
375  }
376  // The WLS integral for a given material
377  // will be inserted in the table according to the
378  // position of the material in the material table.
379 
380  theIntegralTable->insertAt(i,aPhysicsOrderedFreeVector);
381  }
382 }
383 
384 // GetMeanFreePath
385 // ---------------
386 //
388  G4double ,
390 {
391  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
392  const G4Material* aMaterial = aTrack.GetMaterial();
393 
394  G4double thePhotonEnergy = aParticle->GetTotalEnergy();
395 
396  G4MaterialPropertiesTable* aMaterialPropertyTable;
397  G4MaterialPropertyVector* AttenuationLengthVector;
398 
399  G4double AttenuationLength = DBL_MAX;
400 
401  aMaterialPropertyTable = aMaterial->GetMaterialPropertiesTable();
402 
403  if ( aMaterialPropertyTable ) {
404  AttenuationLengthVector = aMaterialPropertyTable->
405  GetProperty("WLSABSLENGTH");
406  if ( AttenuationLengthVector ){
407  AttenuationLength = AttenuationLengthVector->
408  Value(thePhotonEnergy);
409  }
410  else {
411  // G4cout << "No WLS absorption length specified" << G4endl;
412  }
413  }
414  else {
415  // G4cout << "No WLS absortion length specified" << G4endl;
416  }
417 
418  return AttenuationLength;
419 }
420 
421 #include "G4AutoLock.hh"
422 namespace {
423  G4Mutex wlsCmdHandlingMutex = G4MUTEX_INITIALIZER;
424 }
425 
427 {
428  G4AutoLock l(&wlsCmdHandlingMutex);
429  if (name == "delta")
430  {
433  new G4WLSTimeGeneratorProfileDelta("delta");
434  }
435  else if (name == "exponential")
436  {
439  new G4WLSTimeGeneratorProfileExponential("exponential");
440  }
441  else
442  {
443  G4Exception("G4OpWLS::UseTimeProfile", "em0202",
445  "generator does not exist");
446  }
447 }
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
G4int GetNumberOfSecondaries() const
void BuildThePhysicsTable()
Definition: G4OpWLS.cc:295
G4MaterialPropertyVector * GetProperty(const char *key)
G4int verboseLevel
Definition: G4VProcess.hh:368
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4double GetTotalEnergy() const
const G4DynamicParticle * GetDynamicParticle() const
G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep)
Definition: G4OpWLS.cc:109
size_t GetIndex() const
Definition: G4Material.hh:260
G4String name
Definition: TRTMaterials.hh:40
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:564
std::vector< G4Material * > G4MaterialTable
static G4VWLSTimeGeneratorProfile * WLSTimeGeneratorProfile
Definition: G4OpWLS.hh:142
void SetTouchableHandle(const G4TouchableHandle &apValue)
size_t GetVectorLength() const
int G4int
Definition: G4Types.hh:78
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:158
static void UseTimeProfile(const G4String name)
Definition: G4OpWLS.cc:426
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
virtual G4double GenerateTime(const G4double time_constant)=0
const G4ThreeVector & GetPosition() const
void BuildPhysicsTable(const G4ParticleDefinition &aParticleType)
Definition: G4OpWLS.cc:100
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:432
~G4OpWLS()
Definition: G4OpWLS.cc:88
Definition: G4Step.hh:76
G4int GetTrackID() const
void SetPolarization(G4double polX, G4double polY, G4double polZ)
G4double Energy(size_t index) const
static size_t GetNumberOfMaterials()
Definition: G4Material.cc:571
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
void SetKineticEnergy(G4double aEnergy)
const G4TouchableHandle & GetTouchableHandle() const
G4double GetEnergy(G4double aValue)
G4Material * GetMaterial() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static G4OpticalPhoton * OpticalPhoton()
virtual void Initialize(const G4Track &)
G4int G4Mutex
Definition: G4Threading.hh:156
G4MaterialPropertiesTable * GetMaterialPropertiesTable() const
Definition: G4Material.hh:250
void SetNumberOfSecondaries(G4int totSecondaries)
void SetParentID(const G4int aValue)
G4StepPoint * GetPostStepPoint() const
G4ParticleChange aParticleChange
Definition: G4VProcess.hh:289
G4bool ConstPropertyExists(const char *key)
G4OpWLS(const G4String &processName="OpWLS", G4ProcessType type=fOptical)
Definition: G4OpWLS.cc:70
void AddSecondary(G4Track *aSecondary)
void insertAt(size_t, G4PhysicsVector *)
#define G4endl
Definition: G4ios.hh:61
G4double GetGlobalTime() const
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
G4ForceCondition
G4PhysicsTable * theIntegralTable
Definition: G4OpWLS.hh:143
G4ThreeVector G4ParticleMomentum
#define DBL_MAX
Definition: templates.hh:83
#define ns
Definition: xmlparse.cc:597
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
G4double GetMeanFreePath(const G4Track &aTrack, G4double, G4ForceCondition *)
Definition: G4OpWLS.cc:387
void clearAndDestroy()
G4ProcessType