Context Navigation

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

Last change on this file since 1340 was 1340, checked in by garnier, 14 years ago
update ti head
File size: 29.1 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.36 2010/11/03 19:23:48 gum Exp $
28	// GEANT4 tag $Name: xrays-V09-03-05 $
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-10-20 Allow the scintillation yield to be a function
40	// of energy deposited by particle type
41	// Thanks to Zach Hartwig (Department of Nuclear
42	// Science and Engineeering - MIT)
43	// 2010-09-22 by Peter Gumplinger
44	// > scintillation rise time included, thanks to
45	// > Martin Goettlich/DESY
46	// 2005-08-17 by Peter Gumplinger
47	// > change variable name MeanNumPhotons -> MeanNumberOfPhotons
48	// 2005-07-28 by Peter Gumplinger
49	// > add G4ProcessType to constructor
50	// 2004-08-05 by Peter Gumplinger
51	// > changed StronglyForced back to Forced in GetMeanLifeTime
52	// 2002-11-21 by Peter Gumplinger
53	// > change to use G4Poisson for small MeanNumberOfPhotons
54	// 2002-11-07 by Peter Gumplinger
55	// > now allow for fast and slow scintillation component
56	// 2002-11-05 by Peter Gumplinger
57	// > now use scintillation constants from G4Material
58	// 2002-05-09 by Peter Gumplinger
59	// > use only the PostStepPoint location for the origin of
60	// scintillation photons when energy is lost to the medium
61	// by a neutral particle
62	// 2000-09-18 by Peter Gumplinger
63	// > change: aSecondaryPosition=x0+rand*aStep.GetDeltaPosition();
64	// aSecondaryTrack->SetTouchable(0);
65	// 2001-09-17, migration of Materials to pure STL (mma)
66	// 2003-06-03, V.Ivanchenko fix compilation warnings
67	//
68	// mail: gum@triumf.ca
69	//
70	////////////////////////////////////////////////////////////////////////
71
72	#include "G4ios.hh"
73	#include "G4ParticleTypes.hh"
74	#include "G4EmProcessSubType.hh"
75
76	#include "G4Scintillation.hh"
77
78	/////////////////////////
79	// Class Implementation
80	/////////////////////////
81
82	//////////////
83	// Operators
84	//////////////
85
86	// G4Scintillation::operator=(const G4Scintillation &right)
87	// {
88	// }
89
90	/////////////////
91	// Constructors
92	/////////////////
93
94	G4Scintillation::G4Scintillation(const G4String& processName,
95	G4ProcessType type)
96	: G4VRestDiscreteProcess(processName, type)
97	{
98	SetProcessSubType(fScintillation);
99
100	fTrackSecondariesFirst = false;
101	fFiniteRiseTime = false;
102
103	YieldFactor = 1.0;
104	ExcitationRatio = 1.0;
105
106	scintillationByParticleType = false;
107
108	theFastIntegralTable = NULL;
109	theSlowIntegralTable = NULL;
110
111	if (verboseLevel>0) {
112	G4cout << GetProcessName() << " is created " << G4endl;
113	}
114
115	BuildThePhysicsTable();
116
117	emSaturation = NULL;
118	}
119
120	////////////////
121	// Destructors
122	////////////////
123
124	G4Scintillation::~G4Scintillation()
125	{
126	if (theFastIntegralTable != NULL) {
127	theFastIntegralTable->clearAndDestroy();
128	delete theFastIntegralTable;
129	}
130	if (theSlowIntegralTable != NULL) {
131	theSlowIntegralTable->clearAndDestroy();
132	delete theSlowIntegralTable;
133	}
134	}
135
136	////////////
137	// Methods
138	////////////
139
140	// AtRestDoIt
141	// ----------
142	//
143	G4VParticleChange*
144	G4Scintillation::AtRestDoIt(const G4Track& aTrack, const G4Step& aStep)
145
146	// This routine simply calls the equivalent PostStepDoIt since all the
147	// necessary information resides in aStep.GetTotalEnergyDeposit()
148
149	{
150	return G4Scintillation::PostStepDoIt(aTrack, aStep);
151	}
152
153	// PostStepDoIt
154	// -------------
155	//
156	G4VParticleChange*
157	G4Scintillation::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
158
159	// This routine is called for each tracking step of a charged particle
160	// in a scintillator. A Poisson/Gauss-distributed number of photons is
161	// generated according to the scintillation yield formula, distributed
162	// evenly along the track segment and uniformly into 4pi.
163
164	{
165	aParticleChange.Initialize(aTrack);
166
167	const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
168	const G4Material* aMaterial = aTrack.GetMaterial();
169
170	G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
171	G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
172
173	G4ThreeVector x0 = pPreStepPoint->GetPosition();
174	G4ThreeVector p0 = aStep.GetDeltaPosition().unit();
175	G4double t0 = pPreStepPoint->GetGlobalTime();
176
177	G4double TotalEnergyDeposit = aStep.GetTotalEnergyDeposit();
178
179	G4MaterialPropertiesTable* aMaterialPropertiesTable =
180	aMaterial->GetMaterialPropertiesTable();
181	if (!aMaterialPropertiesTable)
182	return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
183
184	const G4MaterialPropertyVector* Fast_Intensity =
185	aMaterialPropertiesTable->GetProperty("FASTCOMPONENT");
186	const G4MaterialPropertyVector* Slow_Intensity =
187	aMaterialPropertiesTable->GetProperty("SLOWCOMPONENT");
188
189	if (!Fast_Intensity && !Slow_Intensity )
190	return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
191
192	G4int nscnt = 1;
193	if (Fast_Intensity && Slow_Intensity) nscnt = 2;
194
195	G4double ScintillationYield = 0.;
196
197	if (scintillationByParticleType) {
198	// The scintillation response is a function of the energy
199	// deposited by particle types.
200
201	// Get the definition of the current particle
202	G4ParticleDefinition *pDef = aParticle->GetDefinition();
203	const G4MaterialPropertyVector *Scint_Yield_Vector = NULL;
204
205	// Obtain the G4MaterialPropertyVectory containing the
206	// scintillation light yield as a function of the deposited
207	// energy for the current particle type
208
209	// Protons
210	if(pDef==G4Proton::ProtonDefinition())
211	Scint_Yield_Vector = aMaterialPropertiesTable->
212	GetProperty("PROTONSCINTILLATIONYIELD");
213
214	// Deuterons
215	else if(pDef==G4Deuteron::DeuteronDefinition())
216	Scint_Yield_Vector = aMaterialPropertiesTable->
217	GetProperty("DEUTERONSCINTILLATIONYIELD");
218
219	// Tritons
220	else if(pDef==G4Triton::TritonDefinition())
221	Scint_Yield_Vector = aMaterialPropertiesTable->
222	GetProperty("TRITONSCINTILLATIONYIELD");
223
224	// Alphas
225	else if(pDef==G4Alpha::AlphaDefinition())
226	Scint_Yield_Vector = aMaterialPropertiesTable->
227	GetProperty("ALPHASCINTILLATIONYIELD");
228
229	// Ions (particles derived from G4VIon and G4Ions)
230	// and recoil ions below tracking cut from neutrons after hElastic
231	else if(pDef->GetParticleType()== "nucleus" \|\|
232	pDef==G4Neutron::NeutronDefinition())
233	Scint_Yield_Vector = aMaterialPropertiesTable->
234	GetProperty("IONSCINTILLATIONYIELD");
235
236	// Electrons (must also account for shell-binding energy
237	// attributed to gamma from standard PhotoElectricEffect)
238	else if(pDef==G4Electron::ElectronDefinition() \|\|
239	pDef==G4Gamma::GammaDefinition())
240	Scint_Yield_Vector = aMaterialPropertiesTable->
241	GetProperty("ELECTRONSCINTILLATIONYIELD");
242
243	// Default for particles not enumerated/listed above
244	else
245	Scint_Yield_Vector = aMaterialPropertiesTable->
246	GetProperty("ELECTRONSCINTILLATIONYIELD");
247
248	// If the user has not specified yields for (p,d,t,a,carbon)
249	// then these unspecified particles will default to the
250	// electron's scintillation yield
251	if(!Scint_Yield_Vector){
252	Scint_Yield_Vector = aMaterialPropertiesTable->
253	GetProperty("ELECTRONSCINTILLATIONYIELD");
254	}
255
256	// Throw an exception if no scintillation yield is found
257	if (!Scint_Yield_Vector) {
258	G4cerr << "\nG4Scintillation::PostStepDoIt(): "
259	<< "Request for scintillation yield for energy deposit and particle type without correct entry in MaterialPropertiesTable\n"
260	<< "ScintillationByParticleType requires at minimum that ELECTRONSCINTILLATIONYIELD is set by the user\n"
261	<< G4endl;
262	G4Exception("G4Scintillation::PostStepDoIt",
263	"No correct entry in MaterialPropertiesTable",
264	FatalException,"Missing MaterialPropertiesTable entry.");
265	}
266
267	if (verboseLevel>1) {
268	G4cout << "\n"
269	<< "Particle = " << pDef->GetParticleName() << "\n"
270	<< "Energy Dep. = " << TotalEnergyDeposit/MeV << "\n"
271	<< "Yield = "
272	<< Scint_Yield_Vector->GetProperty(TotalEnergyDeposit)
273	<< "\n" << G4endl;
274	}
275
276	// Obtain the scintillation yield using the total energy
277	// deposited by the particle in this step.
278
279	// Units: [# scintillation photons]
280	ScintillationYield = Scint_Yield_Vector->
281	GetProperty(TotalEnergyDeposit);
282	} else {
283	// The default linear scintillation process
284	G4double ScintillationYield = aMaterialPropertiesTable->
285	GetConstProperty("SCINTILLATIONYIELD");
286
287	// Units: [# scintillation photons / MeV]
288	ScintillationYield *= YieldFactor;
289	}
290
291	G4double ResolutionScale = aMaterialPropertiesTable->
292	GetConstProperty("RESOLUTIONSCALE");
293
294	// Birks law saturation:
295
296	G4double constBirks = 0.0;
297
298	constBirks = aMaterial->GetIonisation()->GetBirksConstant();
299
300	G4double MeanNumberOfPhotons;
301
302	// Birk's correction via emSaturation and specifying scintillation by
303	// by particle type are physically mutually exclusive
304
305	if (scintillationByParticleType)
306	MeanNumberOfPhotons = ScintillationYield;
307	else if (emSaturation)
308	MeanNumberOfPhotons = ScintillationYield*
309	(emSaturation->VisibleEnergyDeposition(&aStep));
310	else
311	MeanNumberOfPhotons = ScintillationYield*TotalEnergyDeposit;
312
313	G4int NumPhotons;
314
315	if (MeanNumberOfPhotons > 10.)
316	{
317	G4double sigma = ResolutionScale * std::sqrt(MeanNumberOfPhotons);
318	NumPhotons = G4int(G4RandGauss::shoot(MeanNumberOfPhotons,sigma)+0.5);
319	}
320	else
321	{
322	NumPhotons = G4int(G4Poisson(MeanNumberOfPhotons));
323	}
324
325	if (NumPhotons <= 0)
326	{
327	// return unchanged particle and no secondaries
328
329	aParticleChange.SetNumberOfSecondaries(0);
330
331	return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
332	}
333
334	////////////////////////////////////////////////////////////////
335
336	aParticleChange.SetNumberOfSecondaries(NumPhotons);
337
338	if (fTrackSecondariesFirst) {
339	if (aTrack.GetTrackStatus() == fAlive )
340	aParticleChange.ProposeTrackStatus(fSuspend);
341	}
342
343	////////////////////////////////////////////////////////////////
344
345	G4int materialIndex = aMaterial->GetIndex();
346
347	// Retrieve the Scintillation Integral for this material
348	// new G4PhysicsOrderedFreeVector allocated to hold CII's
349
350	G4int Num = NumPhotons;
351
352	for (G4int scnt = 1; scnt <= nscnt; scnt++) {
353
354	G4double ScintillationTime = 0.*ns;
355	G4double ScintillationRiseTime = 0.*ns;
356	G4PhysicsOrderedFreeVector* ScintillationIntegral = NULL;
357
358	if (scnt == 1) {
359	if (nscnt == 1) {
360	if(Fast_Intensity){
361	ScintillationTime = aMaterialPropertiesTable->
362	GetConstProperty("FASTTIMECONSTANT");
363	if (fFiniteRiseTime) {
364	ScintillationRiseTime = aMaterialPropertiesTable->
365	GetConstProperty("FASTSCINTILLATIONRISETIME");
366	}
367	ScintillationIntegral =
368	(G4PhysicsOrderedFreeVector)((theFastIntegralTable)(materialIndex));
369	}
370	if(Slow_Intensity){
371	ScintillationTime = aMaterialPropertiesTable->
372	GetConstProperty("SLOWTIMECONSTANT");
373	if (fFiniteRiseTime) {
374	ScintillationRiseTime = aMaterialPropertiesTable->
375	GetConstProperty("SLOWSCINTILLATIONRISETIME");
376	}
377	ScintillationIntegral =
378	(G4PhysicsOrderedFreeVector)((theSlowIntegralTable)(materialIndex));
379	}
380	}
381	else {
382	G4double YieldRatio = aMaterialPropertiesTable->
383	GetConstProperty("YIELDRATIO");
384	if ( ExcitationRatio == 1.0 ) {
385	Num = G4int (std::min(YieldRatio,1.0) * NumPhotons);
386	}
387	else {
388	Num = G4int (std::min(ExcitationRatio,1.0) * NumPhotons);
389	}
390	ScintillationTime = aMaterialPropertiesTable->
391	GetConstProperty("FASTTIMECONSTANT");
392	if (fFiniteRiseTime) {
393	ScintillationRiseTime = aMaterialPropertiesTable->
394	GetConstProperty("FASTSCINTILLATIONRISETIME");
395	}
396	ScintillationIntegral =
397	(G4PhysicsOrderedFreeVector)((theFastIntegralTable)(materialIndex));
398	}
399	}
400	else {
401	Num = NumPhotons - Num;
402	ScintillationTime = aMaterialPropertiesTable->
403	GetConstProperty("SLOWTIMECONSTANT");
404	if (fFiniteRiseTime) {
405	ScintillationRiseTime = aMaterialPropertiesTable->
406	GetConstProperty("SLOWSCINTILLATIONRISETIME");
407	}
408	ScintillationIntegral =
409	(G4PhysicsOrderedFreeVector)((theSlowIntegralTable)(materialIndex));
410	}
411
412	if (!ScintillationIntegral) continue;
413
414	// Max Scintillation Integral
415
416	G4double CIImax = ScintillationIntegral->GetMaxValue();
417
418	for (G4int i = 0; i < Num; i++) {
419
420	// Determine photon energy
421
422	G4double CIIvalue = G4UniformRand()*CIImax;
423	G4double sampledEnergy =
424	ScintillationIntegral->GetEnergy(CIIvalue);
425
426	if (verboseLevel>1) {
427	G4cout << "sampledEnergy = " << sampledEnergy << G4endl;
428	G4cout << "CIIvalue = " << CIIvalue << G4endl;
429	}
430
431	// Generate random photon direction
432
433	G4double cost = 1. - 2.*G4UniformRand();
434	G4double sint = std::sqrt((1.-cost)*(1.+cost));
435
436	G4double phi = twopi*G4UniformRand();
437	G4double sinp = std::sin(phi);
438	G4double cosp = std::cos(phi);
439
440	G4double px = sint*cosp;
441	G4double py = sint*sinp;
442	G4double pz = cost;
443
444	// Create photon momentum direction vector
445
446	G4ParticleMomentum photonMomentum(px, py, pz);
447
448	// Determine polarization of new photon
449
450	G4double sx = cost*cosp;
451	G4double sy = cost*sinp;
452	G4double sz = -sint;
453
454	G4ThreeVector photonPolarization(sx, sy, sz);
455
456	G4ThreeVector perp = photonMomentum.cross(photonPolarization);
457
458	phi = twopi*G4UniformRand();
459	sinp = std::sin(phi);
460	cosp = std::cos(phi);
461
462	photonPolarization = cosp * photonPolarization + sinp * perp;
463
464	photonPolarization = photonPolarization.unit();
465
466	// Generate a new photon:
467
468	G4DynamicParticle* aScintillationPhoton =
469	new G4DynamicParticle(G4OpticalPhoton::OpticalPhoton(),
470	photonMomentum);
471	aScintillationPhoton->SetPolarization
472	(photonPolarization.x(),
473	photonPolarization.y(),
474	photonPolarization.z());
475
476	aScintillationPhoton->SetKineticEnergy(sampledEnergy);
477
478	// Generate new G4Track object:
479
480	G4double rand;
481
482	if (aParticle->GetDefinition()->GetPDGCharge() != 0) {
483	rand = G4UniformRand();
484	} else {
485	rand = 1.0;
486	}
487
488	G4double delta = rand * aStep.GetStepLength();
489	G4double deltaTime = delta /
490	((pPreStepPoint->GetVelocity()+
491	pPostStepPoint->GetVelocity())/2.);
492
493	// emission time distribution
494	if (ScintillationRiseTime==0.0) {
495	deltaTime = deltaTime -
496	ScintillationTime * std::log( G4UniformRand() );
497	} else {
498	deltaTime = deltaTime +
499	sample_time(ScintillationRiseTime, ScintillationTime);
500	}
501
502	G4double aSecondaryTime = t0 + deltaTime;
503
504	G4ThreeVector aSecondaryPosition =
505	x0 + rand * aStep.GetDeltaPosition();
506
507	G4Track* aSecondaryTrack =
508	new G4Track(aScintillationPhoton,aSecondaryTime,aSecondaryPosition);
509
510	aSecondaryTrack->SetTouchableHandle(
511	aStep.GetPreStepPoint()->GetTouchableHandle());
512	// aSecondaryTrack->SetTouchableHandle((G4VTouchable*)0);
513
514	aSecondaryTrack->SetParentID(aTrack.GetTrackID());
515
516	aParticleChange.AddSecondary(aSecondaryTrack);
517
518	}
519	}
520
521	if (verboseLevel>0) {
522	G4cout << "\n Exiting from G4Scintillation::DoIt -- NumberOfSecondaries = "
523	<< aParticleChange.GetNumberOfSecondaries() << G4endl;
524	}
525
526	return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
527	}
528
529	// BuildThePhysicsTable for the scintillation process
530	// --------------------------------------------------
531	//
532
533	void G4Scintillation::BuildThePhysicsTable()
534	{
535	if (theFastIntegralTable && theSlowIntegralTable) return;
536
537	const G4MaterialTable* theMaterialTable =
538	G4Material::GetMaterialTable();
539	G4int numOfMaterials = G4Material::GetNumberOfMaterials();
540
541	// create new physics table
542
543	if(!theFastIntegralTable)theFastIntegralTable = new G4PhysicsTable(numOfMaterials);
544	if(!theSlowIntegralTable)theSlowIntegralTable = new G4PhysicsTable(numOfMaterials);
545
546	// loop for materials
547
548	for (G4int i=0 ; i < numOfMaterials; i++)
549	{
550	G4PhysicsOrderedFreeVector* aPhysicsOrderedFreeVector =
551	new G4PhysicsOrderedFreeVector();
552	G4PhysicsOrderedFreeVector* bPhysicsOrderedFreeVector =
553	new G4PhysicsOrderedFreeVector();
554
555	// Retrieve vector of scintillation wavelength intensity for
556	// the material from the material's optical properties table.
557
558	G4Material* aMaterial = (*theMaterialTable)[i];
559
560	G4MaterialPropertiesTable* aMaterialPropertiesTable =
561	aMaterial->GetMaterialPropertiesTable();
562
563	if (aMaterialPropertiesTable) {
564
565	G4MaterialPropertyVector* theFastLightVector =
566	aMaterialPropertiesTable->GetProperty("FASTCOMPONENT");
567
568	if (theFastLightVector) {
569
570	// Retrieve the first intensity point in vector
571	// of (photon energy, intensity) pairs
572
573	theFastLightVector->ResetIterator();
574	++(*theFastLightVector); // advance to 1st entry
575
576	G4double currentIN = theFastLightVector->
577	GetProperty();
578
579	if (currentIN >= 0.0) {
580
581	// Create first (photon energy, Scintillation
582	// Integral pair
583
584	G4double currentPM = theFastLightVector->
585	GetPhotonEnergy();
586
587	G4double currentCII = 0.0;
588
589	aPhysicsOrderedFreeVector->
590	InsertValues(currentPM , currentCII);
591
592	// Set previous values to current ones prior to loop
593
594	G4double prevPM = currentPM;
595	G4double prevCII = currentCII;
596	G4double prevIN = currentIN;
597
598	// loop over all (photon energy, intensity)
599	// pairs stored for this material
600
601	while(++(*theFastLightVector))
602	{
603	currentPM = theFastLightVector->
604	GetPhotonEnergy();
605
606	currentIN = theFastLightVector->
607	GetProperty();
608
609	currentCII = 0.5 * (prevIN + currentIN);
610
611	currentCII = prevCII +
612	(currentPM - prevPM) * currentCII;
613
614	aPhysicsOrderedFreeVector->
615	InsertValues(currentPM, currentCII);
616
617	prevPM = currentPM;
618	prevCII = currentCII;
619	prevIN = currentIN;
620	}
621
622	}
623	}
624
625	G4MaterialPropertyVector* theSlowLightVector =
626	aMaterialPropertiesTable->GetProperty("SLOWCOMPONENT");
627
628	if (theSlowLightVector) {
629
630	// Retrieve the first intensity point in vector
631	// of (photon energy, intensity) pairs
632
633	theSlowLightVector->ResetIterator();
634	++(*theSlowLightVector); // advance to 1st entry
635
636	G4double currentIN = theSlowLightVector->
637	GetProperty();
638
639	if (currentIN >= 0.0) {
640
641	// Create first (photon energy, Scintillation
642	// Integral pair
643
644	G4double currentPM = theSlowLightVector->
645	GetPhotonEnergy();
646
647	G4double currentCII = 0.0;
648
649	bPhysicsOrderedFreeVector->
650	InsertValues(currentPM , currentCII);
651
652	// Set previous values to current ones prior to loop
653
654	G4double prevPM = currentPM;
655	G4double prevCII = currentCII;
656	G4double prevIN = currentIN;
657
658	// loop over all (photon energy, intensity)
659	// pairs stored for this material
660
661	while(++(*theSlowLightVector))
662	{
663	currentPM = theSlowLightVector->
664	GetPhotonEnergy();
665
666	currentIN=theSlowLightVector->
667	GetProperty();
668
669	currentCII = 0.5 * (prevIN + currentIN);
670
671	currentCII = prevCII +
672	(currentPM - prevPM) * currentCII;
673
674	bPhysicsOrderedFreeVector->
675	InsertValues(currentPM, currentCII);
676
677	prevPM = currentPM;
678	prevCII = currentCII;
679	prevIN = currentIN;
680	}
681
682	}
683	}
684	}
685
686	// The scintillation integral(s) for a given material
687	// will be inserted in the table(s) according to the
688	// position of the material in the material table.
689
690	theFastIntegralTable->insertAt(i,aPhysicsOrderedFreeVector);
691	theSlowIntegralTable->insertAt(i,bPhysicsOrderedFreeVector);
692
693	}
694	}
695
696	// Called by the user to set the scintillation yield as a function
697	// of energy deposited by particle type
698
699	void G4Scintillation::SetScintillationByParticleType(const G4bool scintType)
700	{
701	if (emSaturation) {
702	G4Exception("G4Scintillation::SetScintillationByParticleType", "Redefinition",
703	JustWarning, "Birks Saturation is replaced by ScintillationByParticleType!");
704	RemoveSaturation();
705	}
706	scintillationByParticleType = scintType;
707	}
708
709	// GetMeanFreePath
710	// ---------------
711	//
712
713	G4double G4Scintillation::GetMeanFreePath(const G4Track&,
714	G4double ,
715	G4ForceCondition* condition)
716	{
717	*condition = StronglyForced;
718
719	return DBL_MAX;
720
721	}
722
723	// GetMeanLifeTime
724	// ---------------
725	//
726
727	G4double G4Scintillation::GetMeanLifeTime(const G4Track&,
728	G4ForceCondition* condition)
729	{
730	*condition = Forced;
731
732	return DBL_MAX;
733
734	}
735
736	G4double G4Scintillation::sample_time(G4double tau1, G4double tau2)
737	{
738	// tau1: rise time and tau2: decay time
739
740	while(1) {
741	// two random numbers
742	G4double ran1 = G4UniformRand();
743	G4double ran2 = G4UniformRand();
744	//
745	// exponential distribution as envelope function: very efficient
746	//
747	G4double d = (tau1+tau2)/tau2;
748	// make sure the envelope function is
749	// always larger than the bi-exponential
750	G4double t = -1.0tau2std::log(1-ran1);
751	G4double g = d*single_exp(t,tau2);
752	if (ran2 <= bi_exp(t,tau1,tau2)/g) return t;
753	}
754	return -1.0;
755	}

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

Download in other formats: