source: trunk/source/processes/electromagnetic/xrays/src/G4Scintillation.cc@ 960

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

import all except CVS
File size: 18.9 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: G4Scintillation.cc,v 1.26 2006/06/29 19:56:11 gunter Exp $
28// GEANT4 tag $Name: geant4-09-01-patch-02 $
29//
30////////////////////////////////////////////////////////////////////////
31// Scintillation Light Class Implementation
32////////////////////////////////////////////////////////////////////////
33//
34// File: G4Scintillation.cc
35// Description: RestDiscrete Process - Generation of Scintillation Photons
36// Version: 1.0
37// Created: 1998-11-07
38// Author: Peter Gumplinger
39// Updated: 2005-08-17 by Peter Gumplinger
40// > change variable name MeanNumPhotons -> MeanNumberOfPhotons
41// 2005-07-28 by Peter Gumplinger
42// > add G4ProcessType to constructor
43// 2004-08-05 by Peter Gumplinger
44// > changed StronglyForced back to Forced in GetMeanLifeTime
45// 2002-11-21 by Peter Gumplinger
46// > change to use G4Poisson for small MeanNumberOfPhotons
47// 2002-11-07 by Peter Gumplinger
48// > now allow for fast and slow scintillation component
49// 2002-11-05 by Peter Gumplinger
50// > now use scintillation constants from G4Material
51// 2002-05-09 by Peter Gumplinger
52// > use only the PostStepPoint location for the origin of
53// scintillation photons when energy is lost to the medium
54// by a neutral particle
55// 2000-09-18 by Peter Gumplinger
56// > change: aSecondaryPosition=x0+rand*aStep.GetDeltaPosition();
57// aSecondaryTrack->SetTouchable(0);
58// 2001-09-17, migration of Materials to pure STL (mma)
59// 2003-06-03, V.Ivanchenko fix compilation warnings
60//
61// mail: gum@triumf.ca
62//
63////////////////////////////////////////////////////////////////////////
64
65#include "G4ios.hh"
66#include "G4Scintillation.hh"
67
68using namespace std;
69
70/////////////////////////
71// Class Implementation
72/////////////////////////
73
74 //////////////
75 // Operators
76 //////////////
77
78// G4Scintillation::operator=(const G4Scintillation &right)
79// {
80// }
81
82 /////////////////
83 // Constructors
84 /////////////////
85
86G4Scintillation::G4Scintillation(const G4String& processName,
87 G4ProcessType type)
88 : G4VRestDiscreteProcess(processName, type)
89{
90 fTrackSecondariesFirst = false;
91
92 YieldFactor = 1.0;
93 ExcitationRatio = 1.0;
94
95 theFastIntegralTable = NULL;
96 theSlowIntegralTable = NULL;
97
98 if (verboseLevel>0) {
99 G4cout << GetProcessName() << " is created " << G4endl;
100 }
101
102 BuildThePhysicsTable();
103}
104
105 ////////////////
106 // Destructors
107 ////////////////
108
109G4Scintillation::~G4Scintillation()
110{
111 if (theFastIntegralTable != NULL) {
112 theFastIntegralTable->clearAndDestroy();
113 delete theFastIntegralTable;
114 }
115 if (theSlowIntegralTable != NULL) {
116 theSlowIntegralTable->clearAndDestroy();
117 delete theSlowIntegralTable;
118 }
119}
120
121 ////////////
122 // Methods
123 ////////////
124
125// AtRestDoIt
126// ----------
127//
128G4VParticleChange*
129G4Scintillation::AtRestDoIt(const G4Track& aTrack, const G4Step& aStep)
130
131// This routine simply calls the equivalent PostStepDoIt since all the
132// necessary information resides in aStep.GetTotalEnergyDeposit()
133
134{
135 return G4Scintillation::PostStepDoIt(aTrack, aStep);
136}
137
138// PostStepDoIt
139// -------------
140//
141G4VParticleChange*
142G4Scintillation::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
143
144// This routine is called for each tracking step of a charged particle
145// in a scintillator. A Poisson/Gauss-distributed number of photons is
146// generated according to the scintillation yield formula, distributed
147// evenly along the track segment and uniformly into 4pi.
148
149{
150 aParticleChange.Initialize(aTrack);
151
152 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
153 const G4Material* aMaterial = aTrack.GetMaterial();
154
155 G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
156 G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
157
158 G4ThreeVector x0 = pPreStepPoint->GetPosition();
159 G4ThreeVector p0 = aStep.GetDeltaPosition().unit();
160 G4double t0 = pPreStepPoint->GetGlobalTime();
161
162 G4double TotalEnergyDeposit = aStep.GetTotalEnergyDeposit();
163
164 G4MaterialPropertiesTable* aMaterialPropertiesTable =
165 aMaterial->GetMaterialPropertiesTable();
166 if (!aMaterialPropertiesTable)
167 return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
168
169 const G4MaterialPropertyVector* Fast_Intensity =
170 aMaterialPropertiesTable->GetProperty("FASTCOMPONENT");
171 const G4MaterialPropertyVector* Slow_Intensity =
172 aMaterialPropertiesTable->GetProperty("SLOWCOMPONENT");
173
174 if (!Fast_Intensity && !Slow_Intensity )
175 return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
176
177 G4int nscnt = 1;
178 if (Fast_Intensity && Slow_Intensity) nscnt = 2;
179
180 G4double ScintillationYield = aMaterialPropertiesTable->
181 GetConstProperty("SCINTILLATIONYIELD");
182 G4double ResolutionScale = aMaterialPropertiesTable->
183 GetConstProperty("RESOLUTIONSCALE");
184
185 ScintillationYield = YieldFactor * ScintillationYield;
186
187 G4double MeanNumberOfPhotons = ScintillationYield * TotalEnergyDeposit;
188
189 G4int NumPhotons;
190 if (MeanNumberOfPhotons > 10.) {
191 G4double sigma = ResolutionScale * sqrt(MeanNumberOfPhotons);
192 NumPhotons = G4int(G4RandGauss::shoot(MeanNumberOfPhotons,sigma)+0.5);
193 }
194 else {
195 NumPhotons = G4int(G4Poisson(MeanNumberOfPhotons));
196 }
197
198 if (NumPhotons <= 0) {
199
200 // return unchanged particle and no secondaries
201
202 aParticleChange.SetNumberOfSecondaries(0);
203
204 return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
205 }
206
207 ////////////////////////////////////////////////////////////////
208
209 aParticleChange.SetNumberOfSecondaries(NumPhotons);
210
211 if (fTrackSecondariesFirst) {
212 if (aTrack.GetTrackStatus() == fAlive )
213 aParticleChange.ProposeTrackStatus(fSuspend);
214 }
215
216 ////////////////////////////////////////////////////////////////
217
218 G4int materialIndex = aMaterial->GetIndex();
219
220 // Retrieve the Scintillation Integral for this material
221 // new G4PhysicsOrderedFreeVector allocated to hold CII's
222
223 G4int Num = NumPhotons;
224
225 for (G4int scnt = 1; scnt <= nscnt; scnt++) {
226
227 G4double ScintillationTime = 0.*ns;
228 G4PhysicsOrderedFreeVector* ScintillationIntegral = NULL;
229
230 if (scnt == 1) {
231 if (nscnt == 1) {
232 if(Fast_Intensity){
233 ScintillationTime = aMaterialPropertiesTable->
234 GetConstProperty("FASTTIMECONSTANT");
235 ScintillationIntegral =
236 (G4PhysicsOrderedFreeVector*)((*theFastIntegralTable)(materialIndex));
237 }
238 if(Slow_Intensity){
239 ScintillationTime = aMaterialPropertiesTable->
240 GetConstProperty("SLOWTIMECONSTANT");
241 ScintillationIntegral =
242 (G4PhysicsOrderedFreeVector*)((*theSlowIntegralTable)(materialIndex));
243 }
244 }
245 else {
246 G4double YieldRatio = aMaterialPropertiesTable->
247 GetConstProperty("YIELDRATIO");
248 if ( ExcitationRatio == 1.0 ) {
249 Num = G4int (min(YieldRatio,1.0) * NumPhotons);
250 }
251 else {
252 Num = G4int (min(ExcitationRatio,1.0) * NumPhotons);
253 }
254 ScintillationTime = aMaterialPropertiesTable->
255 GetConstProperty("FASTTIMECONSTANT");
256 ScintillationIntegral =
257 (G4PhysicsOrderedFreeVector*)((*theFastIntegralTable)(materialIndex));
258 }
259 }
260 else {
261 Num = NumPhotons - Num;
262 ScintillationTime = aMaterialPropertiesTable->
263 GetConstProperty("SLOWTIMECONSTANT");
264 ScintillationIntegral =
265 (G4PhysicsOrderedFreeVector*)((*theSlowIntegralTable)(materialIndex));
266 }
267
268 if (!ScintillationIntegral) continue;
269
270 // Max Scintillation Integral
271
272 G4double CIImax = ScintillationIntegral->GetMaxValue();
273
274 for (G4int i = 0; i < Num; i++) {
275
276 // Determine photon momentum
277
278 G4double CIIvalue = G4UniformRand()*CIImax;
279 G4double sampledMomentum =
280 ScintillationIntegral->GetEnergy(CIIvalue);
281
282 if (verboseLevel>1) {
283 G4cout << "sampledMomentum = " << sampledMomentum << G4endl;
284 G4cout << "CIIvalue = " << CIIvalue << G4endl;
285 }
286
287 // Generate random photon direction
288
289 G4double cost = 1. - 2.*G4UniformRand();
290 G4double sint = sqrt((1.-cost)*(1.+cost));
291
292 G4double phi = twopi*G4UniformRand();
293 G4double sinp = sin(phi);
294 G4double cosp = cos(phi);
295
296 G4double px = sint*cosp;
297 G4double py = sint*sinp;
298 G4double pz = cost;
299
300 // Create photon momentum direction vector
301
302 G4ParticleMomentum photonMomentum(px, py, pz);
303
304 // Determine polarization of new photon
305
306 G4double sx = cost*cosp;
307 G4double sy = cost*sinp;
308 G4double sz = -sint;
309
310 G4ThreeVector photonPolarization(sx, sy, sz);
311
312 G4ThreeVector perp = photonMomentum.cross(photonPolarization);
313
314 phi = twopi*G4UniformRand();
315 sinp = sin(phi);
316 cosp = cos(phi);
317
318 photonPolarization = cosp * photonPolarization + sinp * perp;
319
320 photonPolarization = photonPolarization.unit();
321
322 // Generate a new photon:
323
324 G4DynamicParticle* aScintillationPhoton =
325 new G4DynamicParticle(G4OpticalPhoton::OpticalPhoton(),
326 photonMomentum);
327 aScintillationPhoton->SetPolarization
328 (photonPolarization.x(),
329 photonPolarization.y(),
330 photonPolarization.z());
331
332 aScintillationPhoton->SetKineticEnergy(sampledMomentum);
333
334 // Generate new G4Track object:
335
336 G4double rand;
337
338 if (aParticle->GetDefinition()->GetPDGCharge() != 0) {
339 rand = G4UniformRand();
340 } else {
341 rand = 1.0;
342 }
343
344 G4double delta = rand * aStep.GetStepLength();
345 G4double deltaTime = delta /
346 ((pPreStepPoint->GetVelocity()+
347 pPostStepPoint->GetVelocity())/2.);
348
349 deltaTime = deltaTime -
350 ScintillationTime * log( G4UniformRand() );
351
352 G4double aSecondaryTime = t0 + deltaTime;
353
354 G4ThreeVector aSecondaryPosition =
355 x0 + rand * aStep.GetDeltaPosition();
356
357 G4Track* aSecondaryTrack =
358 new G4Track(aScintillationPhoton,aSecondaryTime,aSecondaryPosition);
359
360 aSecondaryTrack->SetTouchableHandle((G4VTouchable*)0);
361
362 aSecondaryTrack->SetParentID(aTrack.GetTrackID());
363
364 aParticleChange.AddSecondary(aSecondaryTrack);
365
366 }
367 }
368
369 if (verboseLevel>0) {
370 G4cout << "\n Exiting from G4Scintillation::DoIt -- NumberOfSecondaries = "
371 << aParticleChange.GetNumberOfSecondaries() << G4endl;
372 }
373
374 return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
375}
376
377// BuildThePhysicsTable for the scintillation process
378// --------------------------------------------------
379//
380
381void G4Scintillation::BuildThePhysicsTable()
382{
383 if (theFastIntegralTable && theSlowIntegralTable) return;
384
385 const G4MaterialTable* theMaterialTable =
386 G4Material::GetMaterialTable();
387 G4int numOfMaterials = G4Material::GetNumberOfMaterials();
388
389 // create new physics table
390
391 if(!theFastIntegralTable)theFastIntegralTable = new G4PhysicsTable(numOfMaterials);
392 if(!theSlowIntegralTable)theSlowIntegralTable = new G4PhysicsTable(numOfMaterials);
393
394 // loop for materials
395
396 for (G4int i=0 ; i < numOfMaterials; i++)
397 {
398 G4PhysicsOrderedFreeVector* aPhysicsOrderedFreeVector =
399 new G4PhysicsOrderedFreeVector();
400 G4PhysicsOrderedFreeVector* bPhysicsOrderedFreeVector =
401 new G4PhysicsOrderedFreeVector();
402
403 // Retrieve vector of scintillation wavelength intensity for
404 // the material from the material's optical properties table.
405
406 G4Material* aMaterial = (*theMaterialTable)[i];
407
408 G4MaterialPropertiesTable* aMaterialPropertiesTable =
409 aMaterial->GetMaterialPropertiesTable();
410
411 if (aMaterialPropertiesTable) {
412
413 G4MaterialPropertyVector* theFastLightVector =
414 aMaterialPropertiesTable->GetProperty("FASTCOMPONENT");
415
416 if (theFastLightVector) {
417
418 // Retrieve the first intensity point in vector
419 // of (photon momentum, intensity) pairs
420
421 theFastLightVector->ResetIterator();
422 ++(*theFastLightVector); // advance to 1st entry
423
424 G4double currentIN = theFastLightVector->
425 GetProperty();
426
427 if (currentIN >= 0.0) {
428
429 // Create first (photon momentum, Scintillation
430 // Integral pair
431
432 G4double currentPM = theFastLightVector->
433 GetPhotonMomentum();
434
435 G4double currentCII = 0.0;
436
437 aPhysicsOrderedFreeVector->
438 InsertValues(currentPM , currentCII);
439
440 // Set previous values to current ones prior to loop
441
442 G4double prevPM = currentPM;
443 G4double prevCII = currentCII;
444 G4double prevIN = currentIN;
445
446 // loop over all (photon momentum, intensity)
447 // pairs stored for this material
448
449 while(++(*theFastLightVector))
450 {
451 currentPM = theFastLightVector->
452 GetPhotonMomentum();
453
454 currentIN=theFastLightVector->
455 GetProperty();
456
457 currentCII = 0.5 * (prevIN + currentIN);
458
459 currentCII = prevCII +
460 (currentPM - prevPM) * currentCII;
461
462 aPhysicsOrderedFreeVector->
463 InsertValues(currentPM, currentCII);
464
465 prevPM = currentPM;
466 prevCII = currentCII;
467 prevIN = currentIN;
468 }
469
470 }
471 }
472
473 G4MaterialPropertyVector* theSlowLightVector =
474 aMaterialPropertiesTable->GetProperty("SLOWCOMPONENT");
475
476 if (theSlowLightVector) {
477
478 // Retrieve the first intensity point in vector
479 // of (photon momentum, intensity) pairs
480
481 theSlowLightVector->ResetIterator();
482 ++(*theSlowLightVector); // advance to 1st entry
483
484 G4double currentIN = theSlowLightVector->
485 GetProperty();
486
487 if (currentIN >= 0.0) {
488
489 // Create first (photon momentum, Scintillation
490 // Integral pair
491
492 G4double currentPM = theSlowLightVector->
493 GetPhotonMomentum();
494
495 G4double currentCII = 0.0;
496
497 bPhysicsOrderedFreeVector->
498 InsertValues(currentPM , currentCII);
499
500 // Set previous values to current ones prior to loop
501
502 G4double prevPM = currentPM;
503 G4double prevCII = currentCII;
504 G4double prevIN = currentIN;
505
506 // loop over all (photon momentum, intensity)
507 // pairs stored for this material
508
509 while(++(*theSlowLightVector))
510 {
511 currentPM = theSlowLightVector->
512 GetPhotonMomentum();
513
514 currentIN=theSlowLightVector->
515 GetProperty();
516
517 currentCII = 0.5 * (prevIN + currentIN);
518
519 currentCII = prevCII +
520 (currentPM - prevPM) * currentCII;
521
522 bPhysicsOrderedFreeVector->
523 InsertValues(currentPM, currentCII);
524
525 prevPM = currentPM;
526 prevCII = currentCII;
527 prevIN = currentIN;
528 }
529
530 }
531 }
532 }
533
534 // The scintillation integral(s) for a given material
535 // will be inserted in the table(s) according to the
536 // position of the material in the material table.
537
538 theFastIntegralTable->insertAt(i,aPhysicsOrderedFreeVector);
539 theSlowIntegralTable->insertAt(i,bPhysicsOrderedFreeVector);
540
541 }
542}
543
544// GetMeanFreePath
545// ---------------
546//
547
548G4double G4Scintillation::GetMeanFreePath(const G4Track&,
549 G4double ,
550 G4ForceCondition* condition)
551{
552 *condition = StronglyForced;
553
554 return DBL_MAX;
555
556}
557
558// GetMeanLifeTime
559// ---------------
560//
561
562G4double G4Scintillation::GetMeanLifeTime(const G4Track&,
563 G4ForceCondition* condition)
564{
565 *condition = Forced;
566
567 return DBL_MAX;
568
569}
Note: See TracBrowser for help on using the repository browser.