Changeset 961 for trunk/source/processes/electromagnetic/xrays/src/G4Scintillation.cc

Timestamp:

Apr 6, 2009, 12:21:12 PM (15 years ago)

Author:

garnier

Message:

update processes

File:

• trunk/source/processes/electromagnetic/xrays/src/G4Scintillation.cc (modified) (16 diffs)

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

@@ -25 +25 @@
 //
 //
-// $Id: G4Scintillation.cc,v 1.26 2006/06/29 19:56:11 gunter Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4Scintillation.cc,v 1.30 2008/10/22 01:19:11 gum Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 ////////////////////////////////////////////////////////////////////////
@@ -64 +64 @@
 
 #include "G4ios.hh"
+#include "G4EmProcessSubType.hh"
+
 #include "G4Scintillation.hh"
 
@@ -88 +90 @@
                   : G4VRestDiscreteProcess(processName, type)
 {
+        SetProcessSubType(fScintillation);
+
         fTrackSecondariesFirst = false;
 
@@ -101 +105 @@
 
         BuildThePhysicsTable();
+
+        emSaturation = NULL;
 }
 
@@ -180 +186 @@
         G4double ScintillationYield = aMaterialPropertiesTable->
                                       GetConstProperty("SCINTILLATIONYIELD");
+        ScintillationYield *= YieldFactor;
+
         G4double ResolutionScale    = aMaterialPropertiesTable->
                                       GetConstProperty("RESOLUTIONSCALE");
 
-        ScintillationYield = YieldFactor * ScintillationYield;
-
-        G4double MeanNumberOfPhotons = ScintillationYield * TotalEnergyDeposit;
+        // Birks law saturation:
+
+        G4double constBirks = 0.0;
+
+        constBirks = aMaterial->GetIonisation()->GetBirksConstant();
+
+        G4double MeanNumberOfPhotons;
+
+        if (emSaturation) {
+           MeanNumberOfPhotons = ScintillationYield*
+                              (emSaturation->VisibleEnergyDeposition(&aStep));
+        } else {
+           MeanNumberOfPhotons = ScintillationYield*TotalEnergyDeposit;
+        }
 
         G4int NumPhotons;
-        if (MeanNumberOfPhotons > 10.) {
+
+        if (MeanNumberOfPhotons > 10.)
+        {
           G4double sigma = ResolutionScale * sqrt(MeanNumberOfPhotons);
           NumPhotons = G4int(G4RandGauss::shoot(MeanNumberOfPhotons,sigma)+0.5);
         }
-        else {
+        else
+        {
           NumPhotons = G4int(G4Poisson(MeanNumberOfPhotons));
         }
 
-        if (NumPhotons <= 0) {
-
+        if (NumPhotons <= 0)
+        {
            // return unchanged particle and no secondaries 
 
@@ -274 +296 @@
             for (G4int i = 0; i < Num; i++) {
 
-                // Determine photon momentum
+                // Determine photon energy
 
                 G4double CIIvalue = G4UniformRand()*CIImax;
-                G4double sampledMomentum = 
+                G4double sampledEnergy = 
                               ScintillationIntegral->GetEnergy(CIIvalue);
 
                 if (verboseLevel>1) {
-                   G4cout << "sampledMomentum = " << sampledMomentum << G4endl;
+                   G4cout << "sampledEnergy = " << sampledEnergy << G4endl;
                    G4cout << "CIIvalue =        " << CIIvalue << G4endl;
                 }
@@ -330 +352 @@
                                       photonPolarization.z());
 
-                aScintillationPhoton->SetKineticEnergy(sampledMomentum);
+                aScintillationPhoton->SetKineticEnergy(sampledEnergy);
 
                 // Generate new G4Track object:
@@ -358 +380 @@
                 new G4Track(aScintillationPhoton,aSecondaryTime,aSecondaryPosition);
 
-                aSecondaryTrack->SetTouchableHandle((G4VTouchable*)0);
+                aSecondaryTrack->SetTouchableHandle(
+                                 aStep.GetPreStepPoint()->GetTouchableHandle());
+                // aSecondaryTrack->SetTouchableHandle((G4VTouchable*)0);
 
                 aSecondaryTrack->SetParentID(aTrack.GetTrackID());
@@ -417 +441 @@
                 
                       // Retrieve the first intensity point in vector
-                      // of (photon momentum, intensity) pairs 
+                      // of (photon energy, intensity) pairs 
 
                       theFastLightVector->ResetIterator();
@@ -427 +451 @@
                       if (currentIN >= 0.0) {
 
-                         // Create first (photon momentum, Scintillation 
+                         // Create first (photon energy, Scintillation 
                          // Integral pair  
 
                          G4double currentPM = theFastLightVector->
-                                                 GetPhotonMomentum();
+                                                 GetPhotonEnergy();
 
                          G4double currentCII = 0.0;
@@ -444 +468 @@
                          G4double prevIN  = currentIN;
 
-                         // loop over all (photon momentum, intensity)
+                         // loop over all (photon energy, intensity)
                          // pairs stored for this material  
 
@@ -450 +474 @@
                          {
                                 currentPM = theFastLightVector->
-                                                GetPhotonMomentum();
+                                                GetPhotonEnergy();
 
                                 currentIN=theFastLightVector->  
@@ -477 +501 @@
 
                       // Retrieve the first intensity point in vector
-                      // of (photon momentum, intensity) pairs
+                      // of (photon energy, intensity) pairs
 
                       theSlowLightVector->ResetIterator();
@@ -487 +511 @@
                       if (currentIN >= 0.0) {
 
-                         // Create first (photon momentum, Scintillation
+                         // Create first (photon energy, Scintillation
                          // Integral pair
 
                          G4double currentPM = theSlowLightVector->
-                                                 GetPhotonMomentum();
+                                                 GetPhotonEnergy();
 
                          G4double currentCII = 0.0;
@@ -504 +528 @@
                          G4double prevIN  = currentIN;
 
-                         // loop over all (photon momentum, intensity)
+                         // loop over all (photon energy, intensity)
                          // pairs stored for this material
 
@@ -510 +534 @@
                          {
                                 currentPM = theSlowLightVector->
-                                                GetPhotonMomentum();
+                                                GetPhotonEnergy();
 
                                 currentIN=theSlowLightVector->