Context Navigation

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

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: