Context Navigation

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

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

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

Download in other formats: