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

Timestamp:

Nov 5, 2010, 3:45:55 PM (14 years ago)

Author:

garnier

Message:

update ti head

File:

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

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

@@ -25 +25 @@
 //
 //
-// $Id: G4Scintillation.cc,v 1.32 2010/06/16 15:34:15 gcosmo Exp $
-// GEANT4 tag $Name: geant4-09-04-beta-01 $
+// $Id: G4Scintillation.cc,v 1.36 2010/11/03 19:23:48 gum Exp $
+// GEANT4 tag $Name: xrays-V09-03-05 $
 //
 ////////////////////////////////////////////////////////////////////////
@@ -32 +32 @@
 ////////////////////////////////////////////////////////////////////////
 //
-// File:        G4Scintillation.cc 
+// File:        G4Scintillation.cc
 // Description: RestDiscrete Process - Generation of Scintillation Photons
 // Version:     1.0
-// Created:     1998-11-07  
+// Created:     1998-11-07
 // Author:      Peter Gumplinger
-// Updated:     2010-92-22 by Peter Gumplinger
+// Updated:     2010-10-20 Allow the scintillation yield to be a function
+//              of energy deposited by particle type
+//              Thanks to Zach Hartwig (Department of Nuclear
+//              Science and Engineeering - MIT)
+//              2010-09-22 by Peter Gumplinger
 //              > scintillation rise time included, thanks to
 //              > Martin Goettlich/DESY
@@ -59 +63 @@
 //              > change: aSecondaryPosition=x0+rand*aStep.GetDeltaPosition();
 //                        aSecondaryTrack->SetTouchable(0);
-//              2001-09-17, migration of Materials to pure STL (mma) 
+//              2001-09-17, migration of Materials to pure STL (mma)
 //              2003-06-03, V.Ivanchenko fix compilation warnings
 //
@@ -67 +71 @@
 
 #include "G4ios.hh"
+#include "G4ParticleTypes.hh"
 #include "G4EmProcessSubType.hh"
 
@@ -72 +77 @@
 
 /////////////////////////
-// Class Implementation  
+// Class Implementation
 /////////////////////////
 
@@ -93 +98 @@
         SetProcessSubType(fScintillation);
 
-        fTrackSecondariesFirst = false;
+        fTrackSecondariesFirst = false;
         fFiniteRiseTime = false;
 
@@ -99 +104 @@
         ExcitationRatio = 1.0;
 
+        scintillationByParticleType = false;
+
         theFastIntegralTable = NULL;
         theSlowIntegralTable = NULL;
 
-        if (verboseLevel>0) {
+        if (verboseLevel>0) {
            G4cout << GetProcessName() << " is created " << G4endl;
-        }
-
-        BuildThePhysicsTable();
+        }
+
+        BuildThePhysicsTable();
 
         emSaturation = NULL;
@@ -115 +122 @@
         ////////////////
 
-G4Scintillation::~G4Scintillation() 
+G4Scintillation::~G4Scintillation()
 {
         if (theFastIntegralTable != NULL) {
-           theFastIntegralTable->clearAndDestroy();
+           theFastIntegralTable->clearAndDestroy();
            delete theFastIntegralTable;
-        }
+        }
         if (theSlowIntegralTable != NULL) {
            theSlowIntegralTable->clearAndDestroy();
@@ -161 +168 @@
         const G4Material* aMaterial = aTrack.GetMaterial();
 
-        G4StepPoint* pPreStepPoint  = aStep.GetPreStepPoint();
-        G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
-
-        G4ThreeVector x0 = pPreStepPoint->GetPosition();
+        G4StepPoint* pPreStepPoint  = aStep.GetPreStepPoint();
+        G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
+
+        G4ThreeVector x0 = pPreStepPoint->GetPosition();
         G4ThreeVector p0 = aStep.GetDeltaPosition().unit();
-        G4double      t0 = pPreStepPoint->GetGlobalTime();
+        G4double      t0 = pPreStepPoint->GetGlobalTime();
 
         G4double TotalEnergyDeposit = aStep.GetTotalEnergyDeposit();
@@ -175 +182 @@
              return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
 
-        const G4MaterialPropertyVector* Fast_Intensity = 
+        const G4MaterialPropertyVector* Fast_Intensity = 
                 aMaterialPropertiesTable->GetProperty("FASTCOMPONENT"); 
         const G4MaterialPropertyVector* Slow_Intensity =
@@ -186 +193 @@
         if (Fast_Intensity && Slow_Intensity) nscnt = 2;
 
-        G4double ScintillationYield = aMaterialPropertiesTable->
+        G4double ScintillationYield = 0.;
+
+        if (scintillationByParticleType) {
+           // The scintillation response is a function of the energy
+           // deposited by particle types.
+
+           // Get the definition of the current particle
+           G4ParticleDefinition *pDef = aParticle->GetDefinition();
+           const G4MaterialPropertyVector *Scint_Yield_Vector = NULL;
+
+           // Obtain the G4MaterialPropertyVectory containing the
+           // scintillation light yield as a function of the deposited
+           // energy for the current particle type
+
+           // Protons
+           if(pDef==G4Proton::ProtonDefinition()) 
+             Scint_Yield_Vector = aMaterialPropertiesTable->
+               GetProperty("PROTONSCINTILLATIONYIELD");
+
+           // Deuterons
+           else if(pDef==G4Deuteron::DeuteronDefinition())
+             Scint_Yield_Vector = aMaterialPropertiesTable->
+               GetProperty("DEUTERONSCINTILLATIONYIELD");
+
+           // Tritons
+           else if(pDef==G4Triton::TritonDefinition())
+             Scint_Yield_Vector = aMaterialPropertiesTable->
+               GetProperty("TRITONSCINTILLATIONYIELD");
+
+           // Alphas
+           else if(pDef==G4Alpha::AlphaDefinition())
+             Scint_Yield_Vector = aMaterialPropertiesTable->
+               GetProperty("ALPHASCINTILLATIONYIELD");
+          
+           // Ions (particles derived from G4VIon and G4Ions)
+           // and recoil ions below tracking cut from neutrons after hElastic
+           else if(pDef->GetParticleType()== "nucleus" || 
+                   pDef==G4Neutron::NeutronDefinition()) 
+             Scint_Yield_Vector = aMaterialPropertiesTable->
+               GetProperty("IONSCINTILLATIONYIELD");
+
+           // Electrons (must also account for shell-binding energy
+           // attributed to gamma from standard PhotoElectricEffect)
+           else if(pDef==G4Electron::ElectronDefinition() ||
+                   pDef==G4Gamma::GammaDefinition())
+             Scint_Yield_Vector = aMaterialPropertiesTable->
+               GetProperty("ELECTRONSCINTILLATIONYIELD");
+
+           // Default for particles not enumerated/listed above
+           else
+             Scint_Yield_Vector = aMaterialPropertiesTable->
+               GetProperty("ELECTRONSCINTILLATIONYIELD");
+           
+           // If the user has not specified yields for (p,d,t,a,carbon)
+           // then these unspecified particles will default to the 
+           // electron's scintillation yield
+           if(!Scint_Yield_Vector){
+             Scint_Yield_Vector = aMaterialPropertiesTable->
+               GetProperty("ELECTRONSCINTILLATIONYIELD");
+           }
+
+           // Throw an exception if no scintillation yield is found
+           if (!Scint_Yield_Vector) {
+              G4cerr << "\nG4Scintillation::PostStepDoIt(): "
+                     << "Request for scintillation yield for energy deposit and particle type without correct entry in MaterialPropertiesTable\n"
+                     << "ScintillationByParticleType requires at minimum that ELECTRONSCINTILLATIONYIELD is set by the user\n"
+                     << G4endl;
+             G4Exception("G4Scintillation::PostStepDoIt",
+                         "No correct entry in MaterialPropertiesTable",
+                         FatalException,"Missing MaterialPropertiesTable entry.");
+           }
+
+           if (verboseLevel>1) {
+             G4cout << "\n"
+                    << "Particle = " << pDef->GetParticleName() << "\n"
+                    << "Energy Dep. = " << TotalEnergyDeposit/MeV << "\n"
+                    << "Yield = " 
+                    << Scint_Yield_Vector->GetProperty(TotalEnergyDeposit) 
+                    << "\n" << G4endl;
+           }
+
+           // Obtain the scintillation yield using the total energy
+           // deposited by the particle in this step.
+
+           // Units: [# scintillation photons]
+           ScintillationYield = Scint_Yield_Vector->
+                                            GetProperty(TotalEnergyDeposit);
+        } else {
+           // The default linear scintillation process
+           G4double ScintillationYield = aMaterialPropertiesTable->
                                       GetConstProperty("SCINTILLATIONYIELD");
-        ScintillationYield *= YieldFactor;
+
+           // Units: [# scintillation photons / MeV]
+           ScintillationYield *= YieldFactor;
+        }
 
         G4double ResolutionScale    = aMaterialPropertiesTable->
@@ -201 +300 @@
         G4double MeanNumberOfPhotons;
 
-        if (emSaturation) {
+        // Birk's correction via emSaturation and specifying scintillation by
+        // by particle type are physically mutually exclusive
+
+        if (scintillationByParticleType)
+           MeanNumberOfPhotons = ScintillationYield;
+        else if (emSaturation)
            MeanNumberOfPhotons = ScintillationYield*
                               (emSaturation->VisibleEnergyDeposition(&aStep));
-        } else {
+        else
            MeanNumberOfPhotons = ScintillationYield*TotalEnergyDeposit;
-        }
 
         G4int NumPhotons;
@@ -220 +323 @@
         }
 
-        if (NumPhotons <= 0)
+        if (NumPhotons <= 0)
         {
-           // return unchanged particle and no secondaries 
-
-           aParticleChange.SetNumberOfSecondaries(0);
+           // return unchanged particle and no secondaries 
+
+           aParticleChange.SetNumberOfSecondaries(0);
 
            return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
-        }
-
-        ////////////////////////////////////////////////////////////////
-
-        aParticleChange.SetNumberOfSecondaries(NumPhotons);
-
-        if (fTrackSecondariesFirst) {
+        }
+
+        ////////////////////////////////////////////////////////////////
+
+        aParticleChange.SetNumberOfSecondaries(NumPhotons);
+
+        if (fTrackSecondariesFirst) {
            if (aTrack.GetTrackStatus() == fAlive )
-                   aParticleChange.ProposeTrackStatus(fSuspend);
-        }
-        
-        ////////////////////////////////////////////////////////////////
-
-        G4int materialIndex = aMaterial->GetIndex();
-
-        // Retrieve the Scintillation Integral for this material  
-        // new G4PhysicsOrderedFreeVector allocated to hold CII's
+                  aParticleChange.ProposeTrackStatus(fSuspend);
+        }
+
+        ////////////////////////////////////////////////////////////////
+
+        G4int materialIndex = aMaterial->GetIndex();
+
+        // Retrieve the Scintillation Integral for this material  
+        // new G4PhysicsOrderedFreeVector allocated to hold CII's
 
         G4int Num = NumPhotons;
@@ -270 +373 @@
                    if (fFiniteRiseTime) {
                       ScintillationRiseTime = aMaterialPropertiesTable->
-                                  GetConstProperty("SLOWSCINTILLATIONRISETIME");                   }
+                                  GetConstProperty("SLOWSCINTILLATIONRISETIME");
+                   }
                    ScintillationIntegral =
                    (G4PhysicsOrderedFreeVector*)((*theSlowIntegralTable)(materialIndex));
@@ -310 +414 @@
             // Max Scintillation Integral
  
-            G4double CIImax = ScintillationIntegral->GetMaxValue();
-                
-            for (G4int i = 0; i < Num; i++) {
-
-                // Determine photon energy
+            G4double CIImax = ScintillationIntegral->GetMaxValue();
+
+            for (G4int i = 0; i < Num; i++) {
+
+                // Determine photon energy
 
                 G4double CIIvalue = G4UniformRand()*CIImax;
-                G4double sampledEnergy = 
+                G4double sampledEnergy = 
                               ScintillationIntegral->GetEnergy(CIIvalue);
 
-                if (verboseLevel>1) {
+                if (verboseLevel>1) {
                    G4cout << "sampledEnergy = " << sampledEnergy << G4endl;
-                   G4cout << "CIIvalue =        " << CIIvalue << G4endl;
-                }
-
-                // Generate random photon direction
+                   G4cout << "CIIvalue =        " << CIIvalue << G4endl;
+                }
+
+                // Generate random photon direction
 
                 G4double cost = 1. - 2.*G4UniformRand();
                 G4double sint = std::sqrt((1.-cost)*(1.+cost));
 
-                G4double phi = twopi*G4UniformRand();
-                G4double sinp = std::sin(phi);
-                G4double cosp = std::cos(phi);
-
-                G4double px = sint*cosp;
-                G4double py = sint*sinp;
-                G4double pz = cost;
-
-                // Create photon momentum direction vector 
-
-                G4ParticleMomentum photonMomentum(px, py, pz);
-
-                // Determine polarization of new photon 
-
-                G4double sx = cost*cosp;
-                G4double sy = cost*sinp; 
-                G4double sz = -sint;
-
-                G4ThreeVector photonPolarization(sx, sy, sz);
+                G4double phi = twopi*G4UniformRand();
+                G4double sinp = std::sin(phi);
+                G4double cosp = std::cos(phi);
+
+                G4double px = sint*cosp;
+                G4double py = sint*sinp;
+                G4double pz = cost;
+
+                // Create photon momentum direction vector 
+
+                G4ParticleMomentum photonMomentum(px, py, pz);
+
+                // Determine polarization of new photon 
+
+                G4double sx = cost*cosp;
+                G4double sy = cost*sinp; 
+                G4double sz = -sint;
+
+                G4ThreeVector photonPolarization(sx, sy, sz);
 
                 G4ThreeVector perp = photonMomentum.cross(photonPolarization);
 
-                phi = twopi*G4UniformRand();
-                sinp = std::sin(phi);
-                cosp = std::cos(phi);
+                phi = twopi*G4UniformRand();
+                sinp = std::sin(phi);
+                cosp = std::cos(phi);
 
                 photonPolarization = cosp * photonPolarization + sinp * perp;
@@ -364 +468 @@
                 G4DynamicParticle* aScintillationPhoton =
                   new G4DynamicParticle(G4OpticalPhoton::OpticalPhoton(), 
-                                                         photonMomentum);
-                aScintillationPhoton->SetPolarization
-                                     (photonPolarization.x(),
-                                      photonPolarization.y(),
-                                      photonPolarization.z());
-
-                aScintillationPhoton->SetKineticEnergy(sampledEnergy);
+                                                         photonMomentum);
+                aScintillationPhoton->SetPolarization
+                                     (photonPolarization.x(),
+                                      photonPolarization.y(),
+                                      photonPolarization.z());
+
+                aScintillationPhoton->SetKineticEnergy(sampledEnergy);
 
                 // Generate new G4Track object:
@@ -401 +505 @@
                                     x0 + rand * aStep.GetDeltaPosition();
 
-                G4Track* aSecondaryTrack = 
-                new G4Track(aScintillationPhoton,aSecondaryTime,aSecondaryPosition);
+                G4Track* aSecondaryTrack = 
+                new G4Track(aScintillationPhoton,aSecondaryTime,aSecondaryPosition);
 
                 aSecondaryTrack->SetTouchableHandle(
@@ -410 +514 @@
                 aSecondaryTrack->SetParentID(aTrack.GetTrackID());
 
-                aParticleChange.AddSecondary(aSecondaryTrack);
-
-            }
-        }
-
-        if (verboseLevel>0) {
-        G4cout << "\n Exiting from G4Scintillation::DoIt -- NumberOfSecondaries = " 
-             << aParticleChange.GetNumberOfSecondaries() << G4endl;
-        }
-
-        return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
+                aParticleChange.AddSecondary(aSecondaryTrack);
+
+            }
+        }
+
+        if (verboseLevel>0) {
+        G4cout << "\n Exiting from G4Scintillation::DoIt -- NumberOfSecondaries = " 
+               << aParticleChange.GetNumberOfSecondaries() << G4endl;
+        }
+
+        return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
 }
 
@@ -429 +533 @@
 void G4Scintillation::BuildThePhysicsTable()
 {
-        if (theFastIntegralTable && theSlowIntegralTable) return;
-
-        const G4MaterialTable* theMaterialTable = 
+        if (theFastIntegralTable && theSlowIntegralTable) return;
+
+        const G4MaterialTable* theMaterialTable = 
                                G4Material::GetMaterialTable();
-        G4int numOfMaterials = G4Material::GetNumberOfMaterials();
-
-        // create new physics table
+        G4int numOfMaterials = G4Material::GetNumberOfMaterials();
+
+        // create new physics table
         
-        if(!theFastIntegralTable)theFastIntegralTable = new G4PhysicsTable(numOfMaterials);
+        if(!theFastIntegralTable)theFastIntegralTable = new G4PhysicsTable(numOfMaterials);
         if(!theSlowIntegralTable)theSlowIntegralTable = new G4PhysicsTable(numOfMaterials);
 
-        // loop for materials
-
-        for (G4int i=0 ; i < numOfMaterials; i++)
-        {
-                G4PhysicsOrderedFreeVector* aPhysicsOrderedFreeVector =
+        // loop for materials
+
+        for (G4int i=0 ; i < numOfMaterials; i++)
+        {
+                G4PhysicsOrderedFreeVector* aPhysicsOrderedFreeVector =
                                         new G4PhysicsOrderedFreeVector();
                 G4PhysicsOrderedFreeVector* bPhysicsOrderedFreeVector =
                                         new G4PhysicsOrderedFreeVector();
 
-                // Retrieve vector of scintillation wavelength intensity for
+                // Retrieve vector of scintillation wavelength intensity for
                 // the material from the material's optical properties table.
 
-                G4Material* aMaterial = (*theMaterialTable)[i];
-
-                G4MaterialPropertiesTable* aMaterialPropertiesTable =
-                                aMaterial->GetMaterialPropertiesTable();
-
-                if (aMaterialPropertiesTable) {
-
-                   G4MaterialPropertyVector* theFastLightVector = 
-                   aMaterialPropertiesTable->GetProperty("FASTCOMPONENT");
-
-                   if (theFastLightVector) {
-                
-                      // Retrieve the first intensity point in vector
-                      // of (photon energy, intensity) pairs 
-
-                      theFastLightVector->ResetIterator();
-                      ++(*theFastLightVector);  // advance to 1st entry 
-
-                      G4double currentIN = theFastLightVector->
-                                           GetProperty();
-
-                      if (currentIN >= 0.0) {
-
-                         // Create first (photon energy, Scintillation 
+                G4Material* aMaterial = (*theMaterialTable)[i];
+
+                G4MaterialPropertiesTable* aMaterialPropertiesTable =
+                                aMaterial->GetMaterialPropertiesTable();
+
+                if (aMaterialPropertiesTable) {
+
+                   G4MaterialPropertyVector* theFastLightVector = 
+                   aMaterialPropertiesTable->GetProperty("FASTCOMPONENT");
+
+                   if (theFastLightVector) {
+
+                      // Retrieve the first intensity point in vector
+                      // of (photon energy, intensity) pairs 
+
+                      theFastLightVector->ResetIterator();
+                      ++(*theFastLightVector);  // advance to 1st entry 
+
+                      G4double currentIN = theFastLightVector->
+                                               GetProperty();
+
+                      if (currentIN >= 0.0) {
+
+                         // Create first (photon energy, Scintillation 
                          // Integral pair  
 
-                         G4double currentPM = theFastLightVector->
-                                                 GetPhotonEnergy();
-
-                         G4double currentCII = 0.0;
-
-                         aPhysicsOrderedFreeVector->
-                                 InsertValues(currentPM , currentCII);
-
-                         // Set previous values to current ones prior to loop
-
-                         G4double prevPM  = currentPM;
-                         G4double prevCII = currentCII;
-                         G4double prevIN  = currentIN;
-
-                         // loop over all (photon energy, intensity)
-                         // pairs stored for this material  
-
-                         while(++(*theFastLightVector))
-                         {
-                                currentPM = theFastLightVector->
-                                                GetPhotonEnergy();
-
-                                currentIN=theFastLightVector->  
-                                                GetProperty();
-
-                                currentCII = 0.5 * (prevIN + currentIN);
-
-                                currentCII = prevCII +
-                                             (currentPM - prevPM) * currentCII;
-
-                                aPhysicsOrderedFreeVector->
-                                    InsertValues(currentPM, currentCII);
-
-                                prevPM  = currentPM;
-                                prevCII = currentCII;
-                                prevIN  = currentIN;
-                         }
-
-                      }
-                   }
+                         G4double currentPM = theFastLightVector->
+                                                  GetPhotonEnergy();
+
+                         G4double currentCII = 0.0;
+
+                         aPhysicsOrderedFreeVector->
+                                 InsertValues(currentPM , currentCII);
+
+                         // Set previous values to current ones prior to loop
+
+                         G4double prevPM  = currentPM;
+                         G4double prevCII = currentCII;
+                         G4double prevIN  = currentIN;
+
+                         // loop over all (photon energy, intensity)
+                         // pairs stored for this material  
+
+                         while(++(*theFastLightVector))
+                         {
+                                currentPM = theFastLightVector->
+                                                GetPhotonEnergy();
+
+                                currentIN = theFastLightVector->        
+                                                GetProperty();
+
+                                currentCII = 0.5 * (prevIN + currentIN);
+
+                                currentCII = prevCII +
+                                             (currentPM - prevPM) * currentCII;
+
+                                aPhysicsOrderedFreeVector->
+                                    InsertValues(currentPM, currentCII);
+
+                                prevPM  = currentPM;
+                                prevCII = currentCII;
+                                prevIN  = currentIN;
+                         }
+
+                      }
+                   }
 
                    G4MaterialPropertyVector* theSlowLightVector =
@@ -578 +682 @@
                       }
                    }
-                }
-
-        // The scintillation integral(s) for a given material
-        // will be inserted in the table(s) according to the
-        // position of the material in the material table.
-
-        theFastIntegralTable->insertAt(i,aPhysicsOrderedFreeVector);
+                }
+
+        // The scintillation integral(s) for a given material
+        // will be inserted in the table(s) according to the
+        // position of the material in the material table.
+
+        theFastIntegralTable->insertAt(i,aPhysicsOrderedFreeVector);
         theSlowIntegralTable->insertAt(i,bPhysicsOrderedFreeVector);
 
-        }
+        }
+}
+
+// Called by the user to set the scintillation yield as a function
+// of energy deposited by particle type
+
+void G4Scintillation::SetScintillationByParticleType(const G4bool scintType)
+{
+        if (emSaturation) {
+           G4Exception("G4Scintillation::SetScintillationByParticleType", "Redefinition",
+                       JustWarning, "Birks Saturation is replaced by ScintillationByParticleType!");
+           RemoveSaturation();
+        }
+        scintillationByParticleType = scintType;
 }
 
@@ -600 +717 @@
         *condition = StronglyForced;
 
-        return DBL_MAX;
+        return DBL_MAX;
 
 }