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

Last change on this file was 1347, checked in by garnier, 15 years ago

geant4 tag 9.4
File size: 28.3 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.38 2010/12/15 07:39:26 gunter Exp $
28// GEANT4 tag $Name: geant4-09-04-ref-00 $
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
94G4Scintillation::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
124G4Scintillation::~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//
143G4VParticleChange*
144G4Scintillation::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//
156G4VParticleChange*
157G4Scintillation::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 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
533void 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
699void 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
713G4double 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
727G4double G4Scintillation::GetMeanLifeTime(const G4Track&,
728 G4ForceCondition* condition)
729{
730 *condition = Forced;
731
732 return DBL_MAX;
733
734}
735
736G4double 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.0*tau2*std::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.