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

Last change on this file since 846 was 819, checked in by garnier, 16 years ago

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