Context Navigation

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

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

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

Download in other formats: