Index: trunk/source/processes/electromagnetic/xrays/src/G4Scintillation.cc
===================================================================
--- trunk/source/processes/electromagnetic/xrays/src/G4Scintillation.cc	(revision 819)
+++ trunk/source/processes/electromagnetic/xrays/src/G4Scintillation.cc	(revision 961)
@@ -25,6 +25,6 @@
 //
 //
-// $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,4 +64,6 @@
 
 #include "G4ios.hh"
+#include "G4EmProcessSubType.hh"
+
 #include "G4Scintillation.hh"
 
@@ -88,4 +90,6 @@
                   : G4VRestDiscreteProcess(processName, type)
 {
+        SetProcessSubType(fScintillation);
+
 	fTrackSecondariesFirst = false;
 
@@ -101,4 +105,6 @@
 
 	BuildThePhysicsTable();
+
+        emSaturation = NULL;
 }
 
@@ -180,22 +186,38 @@
         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,12 +296,12 @@
 	    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,5 +352,5 @@
 				      photonPolarization.z());
 
-		aScintillationPhoton->SetKineticEnergy(sampledMomentum);
+		aScintillationPhoton->SetKineticEnergy(sampledEnergy);
 
                 // Generate new G4Track object:
@@ -358,5 +380,7 @@
 		new G4Track(aScintillationPhoton,aSecondaryTime,aSecondaryPosition);
 
-                aSecondaryTrack->SetTouchableHandle((G4VTouchable*)0);
+                aSecondaryTrack->SetTouchableHandle(
+                                 aStep.GetPreStepPoint()->GetTouchableHandle());
+                // aSecondaryTrack->SetTouchableHandle((G4VTouchable*)0);
 
                 aSecondaryTrack->SetParentID(aTrack.GetTrackID());
@@ -417,5 +441,5 @@
 		
 		      // Retrieve the first intensity point in vector
-		      // of (photon momentum, intensity) pairs 
+		      // of (photon energy, intensity) pairs 
 
 		      theFastLightVector->ResetIterator();
@@ -427,9 +451,9 @@
 		      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,5 +468,5 @@
                 	 G4double prevIN  = currentIN;
 
-			 // loop over all (photon momentum, intensity)
+			 // loop over all (photon energy, intensity)
 			 // pairs stored for this material  
 
@@ -450,5 +474,5 @@
 			 {
 				currentPM = theFastLightVector->
-						GetPhotonMomentum();
+						GetPhotonEnergy();
 
 				currentIN=theFastLightVector->	
@@ -477,5 +501,5 @@
 
                       // Retrieve the first intensity point in vector
-                      // of (photon momentum, intensity) pairs
+                      // of (photon energy, intensity) pairs
 
                       theSlowLightVector->ResetIterator();
@@ -487,9 +511,9 @@
                       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,5 +528,5 @@
                          G4double prevIN  = currentIN;
 
-                         // loop over all (photon momentum, intensity)
+                         // loop over all (photon energy, intensity)
                          // pairs stored for this material
 
@@ -510,5 +534,5 @@
                          {
                                 currentPM = theSlowLightVector->
-                                                GetPhotonMomentum();
+                                                GetPhotonEnergy();
 
                                 currentIN=theSlowLightVector->
