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

Last change on this file since 1330 was 1315, checked in by garnier, 15 years ago

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 21.4 KB
RevLine 
[819]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//
[1315]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 $
[819]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
[1315]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
[819]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"
[961]69#include "G4EmProcessSubType.hh"
70
[819]71#include "G4Scintillation.hh"
72
73using 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
91G4Scintillation::G4Scintillation(const G4String& processName,
92 G4ProcessType type)
93 : G4VRestDiscreteProcess(processName, type)
94{
[961]95 SetProcessSubType(fScintillation);
96
[819]97 fTrackSecondariesFirst = false;
[1315]98 fFiniteRiseTime = false;
[819]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();
[961]111
112 emSaturation = NULL;
[819]113}
114
115 ////////////////
116 // Destructors
117 ////////////////
118
119G4Scintillation::~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//
138G4VParticleChange*
139G4Scintillation::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//
151G4VParticleChange*
152G4Scintillation::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");
[961]192 ScintillationYield *= YieldFactor;
193
[819]194 G4double ResolutionScale = aMaterialPropertiesTable->
195 GetConstProperty("RESOLUTIONSCALE");
196
[961]197 // Birks law saturation:
[819]198
[961]199 G4double constBirks = 0.0;
[819]200
[961]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
[819]212 G4int NumPhotons;
[961]213
214 if (MeanNumberOfPhotons > 10.)
215 {
[819]216 G4double sigma = ResolutionScale * sqrt(MeanNumberOfPhotons);
217 NumPhotons = G4int(G4RandGauss::shoot(MeanNumberOfPhotons,sigma)+0.5);
218 }
[961]219 else
220 {
[819]221 NumPhotons = G4int(G4Poisson(MeanNumberOfPhotons));
222 }
223
[961]224 if (NumPhotons <= 0)
225 {
[819]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;
[1315]254 G4double ScintillationRiseTime = 0.*ns;
[819]255 G4PhysicsOrderedFreeVector* ScintillationIntegral = NULL;
256
257 if (scnt == 1) {
258 if (nscnt == 1) {
259 if(Fast_Intensity){
260 ScintillationTime = aMaterialPropertiesTable->
261 GetConstProperty("FASTTIMECONSTANT");
[1315]262 if (fFiniteRiseTime) {
263 ScintillationRiseTime = aMaterialPropertiesTable->
264 GetConstProperty("FASTSCINTILLATIONRISETIME");
265 }
[819]266 ScintillationIntegral =
267 (G4PhysicsOrderedFreeVector*)((*theFastIntegralTable)(materialIndex));
268 }
269 if(Slow_Intensity){
270 ScintillationTime = aMaterialPropertiesTable->
271 GetConstProperty("SLOWTIMECONSTANT");
[1315]272 if (fFiniteRiseTime) {
273 ScintillationRiseTime = aMaterialPropertiesTable->
274 GetConstProperty("SLOWSCINTILLATIONRISETIME"); }
[819]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");
[1315]290 if (fFiniteRiseTime) {
291 ScintillationRiseTime = aMaterialPropertiesTable->
292 GetConstProperty("FASTSCINTILLATIONRISETIME");
293 }
[819]294 ScintillationIntegral =
295 (G4PhysicsOrderedFreeVector*)((*theFastIntegralTable)(materialIndex));
296 }
297 }
298 else {
299 Num = NumPhotons - Num;
300 ScintillationTime = aMaterialPropertiesTable->
301 GetConstProperty("SLOWTIMECONSTANT");
[1315]302 if (fFiniteRiseTime) {
303 ScintillationRiseTime = aMaterialPropertiesTable->
304 GetConstProperty("SLOWSCINTILLATIONRISETIME");
305 }
[819]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
[961]318 // Determine photon energy
[819]319
320 G4double CIIvalue = G4UniformRand()*CIImax;
[961]321 G4double sampledEnergy =
[819]322 ScintillationIntegral->GetEnergy(CIIvalue);
323
324 if (verboseLevel>1) {
[961]325 G4cout << "sampledEnergy = " << sampledEnergy << G4endl;
[819]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
[961]374 aScintillationPhoton->SetKineticEnergy(sampledEnergy);
[819]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
[1315]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 }
[819]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
[961]408 aSecondaryTrack->SetTouchableHandle(
409 aStep.GetPreStepPoint()->GetTouchableHandle());
410 // aSecondaryTrack->SetTouchableHandle((G4VTouchable*)0);
[819]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
431void 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
[961]469 // of (photon energy, intensity) pairs
[819]470
471 theFastLightVector->ResetIterator();
472 ++(*theFastLightVector); // advance to 1st entry
473
474 G4double currentIN = theFastLightVector->
475 GetProperty();
476
477 if (currentIN >= 0.0) {
478
[961]479 // Create first (photon energy, Scintillation
[819]480 // Integral pair
481
482 G4double currentPM = theFastLightVector->
[961]483 GetPhotonEnergy();
[819]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
[961]496 // loop over all (photon energy, intensity)
[819]497 // pairs stored for this material
498
499 while(++(*theFastLightVector))
500 {
501 currentPM = theFastLightVector->
[961]502 GetPhotonEnergy();
[819]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
[961]529 // of (photon energy, intensity) pairs
[819]530
531 theSlowLightVector->ResetIterator();
532 ++(*theSlowLightVector); // advance to 1st entry
533
534 G4double currentIN = theSlowLightVector->
535 GetProperty();
536
537 if (currentIN >= 0.0) {
538
[961]539 // Create first (photon energy, Scintillation
[819]540 // Integral pair
541
542 G4double currentPM = theSlowLightVector->
[961]543 GetPhotonEnergy();
[819]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
[961]556 // loop over all (photon energy, intensity)
[819]557 // pairs stored for this material
558
559 while(++(*theSlowLightVector))
560 {
561 currentPM = theSlowLightVector->
[961]562 GetPhotonEnergy();
[819]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
598G4double 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
612G4double G4Scintillation::GetMeanLifeTime(const G4Track&,
613 G4ForceCondition* condition)
614{
615 *condition = Forced;
616
617 return DBL_MAX;
618
619}
[1315]620
621G4double 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.0*tau2*log(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.