/search.css" rel="stylesheet" type="text/css"/> /search.js">
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

Public Member Functions | Protected Attributes | Private Member Functions | Private Attributes
DsG4Scintillation Class Reference

#include <DsG4Scintillation.h>

List of all members.

Public Member Functions

 DsG4Scintillation (const G4String &processName="Scintillation", G4ProcessType type=fElectromagnetic)
 ~DsG4Scintillation ()
G4bool IsApplicable (const G4ParticleDefinition &aParticleType)
G4double GetMeanFreePath (const G4Track &aTrack, G4double, G4ForceCondition *)
G4double GetMeanLifeTime (const G4Track &aTrack, G4ForceCondition *)
G4VParticleChange * PostStepDoIt (const G4Track &aTrack, const G4Step &aStep)
G4VParticleChange * AtRestDoIt (const G4Track &aTrack, const G4Step &aStep)
void SetTrackSecondariesFirst (const G4bool state)
G4bool GetTrackSecondariesFirst () const
void SetScintillationYieldFactor (const G4double yieldfactor)
G4double GetScintillationYieldFactor () const
void SetUseFastMu300nsTrick (const G4bool fastMu300nsTrick)
G4bool GetUseFastMu300nsTrick () const
void SetScintillationExcitationRatio (const G4double excitationratio)
G4double GetScintillationExcitationRatio () const
G4PhysicsTable * GetFastIntegralTable () const
G4PhysicsTable * GetSlowIntegralTable () const
G4PhysicsTable * GetReemissionIntegralTable () const
void DumpPhysicsTable () const
G4double GetPhotonWeight () const
void SetPhotonWeight (G4double weight)
void SetDoReemission (bool tf=true)
bool GetDoReemission ()
void SetDoBothProcess (bool tf=true)
bool GetDoBothProcess ()
G4bool GetApplyPreQE () const
void SetApplyPreQE (G4bool a)
G4double GetPreQE () const
void SetPreQE (G4double a)
void SetBirksConstant1 (double c1)
double GetBirksConstant1 ()
void SetBirksConstant2 (double c2)
double GetBirksConstant2 ()
void SetGammaSlowerTimeConstant (double st)
double GetGammaSlowerTimeConstant ()
void SetGammaSlowerRatio (double sr)
double GetGammaSlowerRatio ()
void SetNeutronSlowerTimeConstant (double st)
double GetNeutronSlowerTimeConstant ()
void SetNeutronSlowerRatio (double sr)
double GetNeutronSlowerRatio ()
void SetAlphaSlowerTimeConstant (double st)
double GetAlphaSlowerTimeConstant ()
void SetAlphaSlowerRatio (double sr)
double GetAlphaSlowerRatio ()
void SetNoOp (bool tf=true)

Protected Attributes

G4PhysicsTable * theSlowIntegralTable
G4PhysicsTable * theFastIntegralTable
G4PhysicsTable * theReemissionIntegralTable
G4bool doReemission
G4bool doBothProcess
double birksConstant1
double birksConstant2
double slowerTimeConstant
double slowerRatio
double gammaSlowerTime
double gammaSlowerRatio
double neutronSlowerTime
double neutronSlowerRatio
double alphaSlowerTime
double alphaSlowerRatio

Private Member Functions

void BuildThePhysicsTable ()

Private Attributes

G4bool fTrackSecondariesFirst
G4double YieldFactor
G4bool FastMu300nsTrick
G4double ExcitationRatio
G4double fPhotonWeight
G4bool fApplyPreQE
G4double fPreQE
bool m_noop

Detailed Description

Definition at line 103 of file DsG4Scintillation.h.


Constructor & Destructor Documentation

DsG4Scintillation::DsG4Scintillation ( const G4String &  processName = "Scintillation",
G4ProcessType  type = fElectromagnetic 
)

Definition at line 99 of file DsG4Scintillation.cc.

    : G4VRestDiscreteProcess(processName, type)
    , doReemission(true)
    , doBothProcess(true)
    , fPhotonWeight(1.0)
    , fApplyPreQE(false)
    , fPreQE(1.)
    , m_noop(false)
{
    SetProcessSubType(fScintillation);
    fTrackSecondariesFirst = false;

    YieldFactor = 1.0;
    ExcitationRatio = 1.0;

    theFastIntegralTable = NULL;
    theSlowIntegralTable = NULL;
    theReemissionIntegralTable = NULL;

    //verboseLevel = 2;
    //G4cout << " DsG4Scintillation set verboseLevel by hand to " << verboseLevel << G4endl;

    if (verboseLevel > 0) {
        G4cout << GetProcessName() << " is created " << G4endl;
    }

    BuildThePhysicsTable();

}
DsG4Scintillation::~DsG4Scintillation ( )

Definition at line 134 of file DsG4Scintillation.cc.

{
    if (theFastIntegralTable != NULL) {
        theFastIntegralTable->clearAndDestroy();
        delete theFastIntegralTable;
    }
    if (theSlowIntegralTable != NULL) {
        theSlowIntegralTable->clearAndDestroy();
        delete theSlowIntegralTable;
    }
    if (theReemissionIntegralTable != NULL) {
        theReemissionIntegralTable->clearAndDestroy();
        delete theReemissionIntegralTable;
    }
}

Member Function Documentation

G4bool DsG4Scintillation::IsApplicable ( const G4ParticleDefinition &  aParticleType) [inline]

Definition at line 301 of file DsG4Scintillation.h.

{
        if (aParticleType.GetParticleName() == "opticalphoton"){
           return true;
        } else {
           return true;
        }
}
G4double DsG4Scintillation::GetMeanFreePath ( const G4Track &  aTrack,
G4double  ,
G4ForceCondition *  condition 
)

Definition at line 977 of file DsG4Scintillation.cc.

{
    *condition = StronglyForced;

    return DBL_MAX;

}
G4double DsG4Scintillation::GetMeanLifeTime ( const G4Track &  aTrack,
G4ForceCondition *  condition 
)

Definition at line 991 of file DsG4Scintillation.cc.

{
    *condition = Forced;

    return DBL_MAX;

}
G4VParticleChange * DsG4Scintillation::PostStepDoIt ( const G4Track &  aTrack,
const G4Step &  aStep 
)

Definition at line 171 of file DsG4Scintillation.cc.

{
    aParticleChange.Initialize(aTrack);

    if (m_noop) {               // do nothing, bail
        aParticleChange.SetNumberOfSecondaries(0);
        return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
    }


    G4String pname="";
    G4ThreeVector vertpos;
    G4double vertenergy=0.0;
    G4double reem_d=0.0;
    G4bool flagReemission= false;
    DsPhotonTrackInfo* reemittedTI=0;
    if (aTrack.GetDefinition() == G4OpticalPhoton::OpticalPhoton()) {
        G4Track *track=aStep.GetTrack();
        G4CompositeTrackInfo* composite=dynamic_cast<G4CompositeTrackInfo*>(track->GetUserInformation());
        reemittedTI = composite?dynamic_cast<DsPhotonTrackInfo*>( composite->GetPhotonTrackInfo() ):0;
        
        const G4VProcess* process = track->GetCreatorProcess();
        if(process) pname = process->GetProcessName();

        if (verboseLevel > 0) { 
          G4cout<<"Optical photon. Process name is " << pname<<G4endl;
        } 
        if(doBothProcess) {
            flagReemission= doReemission
                && aTrack.GetTrackStatus() == fStopAndKill
                && aStep.GetPostStepPoint()->GetStepStatus() != fGeomBoundary;     
        }
        else{
            flagReemission= doReemission
                && aTrack.GetTrackStatus() == fStopAndKill
                && aStep.GetPostStepPoint()->GetStepStatus() != fGeomBoundary
                && pname=="Cerenkov";
        }
        if(verboseLevel > 0) {
            G4cout<<"flag of Reemission is "<<flagReemission<<"!!"<<G4endl;
        }
        if (!flagReemission) {
            return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
        }
    }

    G4double TotalEnergyDeposit = aStep.GetTotalEnergyDeposit();
    if (verboseLevel > 0 ) { 
      G4cout << " TotalEnergyDeposit " << TotalEnergyDeposit 
             << " material " << aTrack.GetMaterial()->GetName() << G4endl;
    }
    if (TotalEnergyDeposit <= 0.0 && !flagReemission) {
        return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
    }

    const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
    const G4String aParticleName = aParticle->GetDefinition()->GetParticleName();
    const G4Material* aMaterial = aTrack.GetMaterial();

    G4MaterialPropertiesTable* aMaterialPropertiesTable =
        aMaterial->GetMaterialPropertiesTable();
    if (!aMaterialPropertiesTable)
        return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);

    G4String FastTimeConstant = "FASTTIMECONSTANT";
    G4String SlowTimeConstant = "SLOWTIMECONSTANT";
    G4String strYieldRatio = "YIELDRATIO";

    
    if (aParticleName == "opticalphoton") {
      FastTimeConstant = "ReemissionFASTTIMECONSTANT";
      SlowTimeConstant = "ReemissionSLOWTIMECONSTANT";
      strYieldRatio = "ReemissionYIELDRATIO";
    }
    else if(aParticleName == "gamma" || aParticleName == "e+" || aParticleName == "e-") {
      FastTimeConstant = "GammaFASTTIMECONSTANT";
      SlowTimeConstant = "GammaSLOWTIMECONSTANT";
      strYieldRatio = "GammaYIELDRATIO";
      slowerTimeConstant = gammaSlowerTime;
      slowerRatio = gammaSlowerRatio;
    }
    else if(aParticleName == "alpha") {
      FastTimeConstant = "AlphaFASTTIMECONSTANT";
      SlowTimeConstant = "AlphaSLOWTIMECONSTANT";
      strYieldRatio = "AlphaYIELDRATIO";
      slowerTimeConstant = alphaSlowerTime;
      slowerRatio = alphaSlowerRatio;
    }
    else {
      FastTimeConstant = "NeutronFASTTIMECONSTANT";
      SlowTimeConstant = "NeutronSLOWTIMECONSTANT";
      strYieldRatio = "NeutronYIELDRATIO";
      slowerTimeConstant = neutronSlowerTime;
      slowerRatio = neutronSlowerRatio;
    }

    const G4MaterialPropertyVector* Fast_Intensity = 
        aMaterialPropertiesTable->GetProperty("FASTCOMPONENT"); 
    const G4MaterialPropertyVector* Slow_Intensity =
        aMaterialPropertiesTable->GetProperty("SLOWCOMPONENT");
    const G4MaterialPropertyVector* Reemission_Prob =
        aMaterialPropertiesTable->GetProperty("REEMISSIONPROB");
    if (verboseLevel > 0 ) {
      G4cout << " MaterialPropertyVectors: Fast_Intensity " << Fast_Intensity 
             << " Slow_Intensity " << Slow_Intensity << " Reemission_Prob " << Reemission_Prob << G4endl;
    }
    if (!Fast_Intensity && !Slow_Intensity )
        return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);

    G4int nscnt = 1;
    if (Fast_Intensity && Slow_Intensity) nscnt = 2;
    if ( verboseLevel > 0) {
      G4cout << " Fast_Intensity " << Fast_Intensity << " Slow_Intensity " << Slow_Intensity << " nscnt " << nscnt << G4endl;
    }
    G4StepPoint* pPreStepPoint  = aStep.GetPreStepPoint();
    G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();

    G4ThreeVector x0 = pPreStepPoint->GetPosition();
    G4ThreeVector p0 = aStep.GetDeltaPosition().unit();
    G4double      t0 = pPreStepPoint->GetGlobalTime();

    //Replace NumPhotons by NumTracks
    G4int NumTracks=0;
    G4double weight=1.0;
    if (flagReemission) {   
        if(verboseLevel > 0){   
            G4cout<<"the process name is "<<pname<<"!!"<<G4endl;}
        
        if ( Reemission_Prob == 0)
            return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
        G4double p_reemission=
            Reemission_Prob->GetProperty(aTrack.GetKineticEnergy());
        if (G4UniformRand() >= p_reemission)
            return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
        NumTracks= 1;
        weight= aTrack.GetWeight();
        if (verboseLevel > 0 ) {
            G4cout << " flagReemission " << flagReemission << " weight " << weight << G4endl;}
    }
    else {
        // J.B.Birks. The theory and practice of Scintillation Counting. 
        // Pergamon Press, 1964.      
        // For particles with energy much smaller than minimum ionization 
        // energy, the scintillation response is non-linear because of quenching  
        // effect. The light output is reduced by a parametric factor: 
        // 1/(1 + birk1*delta + birk2* delta^2). 
        // Delta is the energy loss per unit mass thickness. birk1 and birk2 
        // were measured for several organic scintillators.         
        // Here we use birk1 = 0.0125*g/cm2/MeV and ignore birk2.               
        // R.L.Craun and D.L.Smith. Nucl. Inst. and Meth., 80:239-244, 1970.   
        // Liang Zhan  01/27/2006 
        // /////////////////////////////////////////////////////////////////////
        
        
        G4double ScintillationYield = 0;
        {// Yield.  Material must have this or we lack raisins dayetras
            const G4MaterialPropertyVector* ptable =
                aMaterialPropertiesTable->GetProperty("SCINTILLATIONYIELD");
            if (!ptable) {
                G4cout << "ConstProperty: failed to get SCINTILLATIONYIELD"
                       << G4endl;
                return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
            }
            ScintillationYield = ptable->GetProperty(0);
        }

        G4double ResolutionScale    = 1;
        {// Resolution Scale
            const G4MaterialPropertyVector* ptable =
                aMaterialPropertiesTable->GetProperty("RESOLUTIONSCALE");
            if (ptable)
                ResolutionScale = ptable->GetProperty(0);
        }

        G4double dE = TotalEnergyDeposit;
        G4double dx = aStep.GetStepLength();
        G4double dE_dx = dE/dx;
        if(aTrack.GetDefinition() == G4Gamma::Gamma() && dE > 0)
        { 
          G4LossTableManager* manager = G4LossTableManager::Instance();
          dE_dx = dE/manager->GetRange(G4Electron::Electron(), dE, aTrack.GetMaterialCutsCouple());
          //G4cout<<"gamma dE_dx = "<<dE_dx/(MeV/mm)<<"MeV/mm"<<G4endl;
        }
        
        G4double delta = dE_dx/aMaterial->GetDensity();//get scintillator density 
        //G4double birk1 = 0.0125*g/cm2/MeV;
        G4double birk1 = birksConstant1;
        if(abs(aParticle->GetCharge())>1.5)//for particle charge greater than 1.
            birk1 = 0.57*birk1;
        
        G4double birk2 = 0;
        //birk2 = (0.0031*g/MeV/cm2)*(0.0031*g/MeV/cm2);
        birk2 = birksConstant2;
        
        G4double QuenchedTotalEnergyDeposit 
            = TotalEnergyDeposit/(1+birk1*delta+birk2*delta*delta);

       //Add 300ns trick for muon simuation, by Haoqi Jan 27, 2011  
       if(FastMu300nsTrick)  {
           // cout<<"GlobalTime ="<<aStep.GetTrack()->GetGlobalTime()/ns<<endl;
           if(aStep.GetTrack()->GetGlobalTime()/ns>300) {
               ScintillationYield = YieldFactor * ScintillationYield;
           }
           else{
            ScintillationYield=0.;
           }
        }
        else {    
            ScintillationYield = YieldFactor * ScintillationYield; 
        }

        G4double MeanNumberOfPhotons= ScintillationYield * QuenchedTotalEnergyDeposit;
   
        // Implemented the fast simulation method from GLG4Scint
        // Jianglai 09-05-2006
        
        // randomize number of TRACKS (not photons)
        // this gets statistics right for number of PE after applying
        // boolean random choice to final absorbed track (change from
        // old method of applying binomial random choice to final absorbed
        // track, which did want poissonian number of photons divided
        // as evenly as possible into tracks)
        // Note for weight=1, there's no difference between tracks and photons.
        G4double MeanNumberOfTracks= MeanNumberOfPhotons/fPhotonWeight; 
        if ( fApplyPreQE ) {
            MeanNumberOfTracks*=fPreQE;
        }
        if (MeanNumberOfTracks > 10.) {
            G4double sigma = ResolutionScale * sqrt(MeanNumberOfTracks);
            NumTracks = G4int(G4RandGauss::shoot(MeanNumberOfTracks,sigma)+0.5);
        }
        else {
            NumTracks = G4int(G4Poisson(MeanNumberOfTracks));
        }
        if ( verboseLevel > 0 ) {
          G4cout << " Generated " << NumTracks << " scint photons. mean(scint photons) = " << MeanNumberOfTracks << G4endl;
        }
    }
    weight*=fPhotonWeight;
    if ( verboseLevel > 0 ) {
      G4cout << " set scint photon weight to " << weight << " after multiplying original weight by fPhotonWeight " << fPhotonWeight 
             << " NumTracks = " << NumTracks
             << G4endl;
    }
    // G4cerr<<"Scint weight is "<<weight<<G4endl;
    if (NumTracks <= 0) {
        // return unchanged particle and no secondaries 
        aParticleChange.SetNumberOfSecondaries(0);
        return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
    }


    aParticleChange.SetNumberOfSecondaries(NumTracks);

    if (fTrackSecondariesFirst) {
        if (!flagReemission) 
            if (aTrack.GetTrackStatus() == fAlive )
                aParticleChange.ProposeTrackStatus(fSuspend);
    }
        

    G4int materialIndex = aMaterial->GetIndex();

    G4PhysicsOrderedFreeVector* ReemissionIntegral = NULL;
    ReemissionIntegral =
        (G4PhysicsOrderedFreeVector*)((*theReemissionIntegralTable)(materialIndex));

    // Retrieve the Scintillation Integral for this material  
    // new G4PhysicsOrderedFreeVector allocated to hold CII's

    G4int Num = NumTracks; //# tracks is now the loop control
        
    G4double fastTimeConstant = 0.0;
    { // Fast Time Constant
        const G4MaterialPropertyVector* ptable =
        aMaterialPropertiesTable->GetProperty(FastTimeConstant.c_str());
        if (verboseLevel > 0) {
          G4cout << " MaterialPropertyVector table " << ptable << " for FASTTIMECONSTANT"<<G4endl;
        }
        if (!ptable) ptable = aMaterialPropertiesTable->GetProperty("FASTTIMECONSTANT");
        if (ptable) {
            fastTimeConstant = ptable->GetProperty(0);
          if (verboseLevel > 0) { 
            G4cout << " dump fast time constant table " << G4endl;
            const_cast <G4MaterialPropertyVector*>(ptable)->DumpVector();
          }
        }
    }

    G4double slowTimeConstant = 0.0;
    { // Slow Time Constant
        const G4MaterialPropertyVector* ptable =
        aMaterialPropertiesTable->GetProperty(SlowTimeConstant.c_str());
        if (verboseLevel > 0) {
          G4cout << " MaterialPropertyVector table " << ptable << " for SLOWTIMECONSTANT"<<G4endl;
        }
        if(!ptable) ptable = aMaterialPropertiesTable->GetProperty("SLOWTIMECONSTANT");
        if (ptable){
          slowTimeConstant = ptable->GetProperty(0);
          if (verboseLevel > 0) { 
            G4cout << " dump slow time constant table " << G4endl;
            const_cast <G4MaterialPropertyVector*>(ptable)->DumpVector();
          }
        }
    }

    G4double YieldRatio = 0.0;
    { // Slow Time Constant
        const G4MaterialPropertyVector* ptable =
            aMaterialPropertiesTable->GetProperty(strYieldRatio.c_str());
        if(!ptable) ptable = aMaterialPropertiesTable->GetProperty("YIELDRATIO");
        if (ptable)
            YieldRatio = ptable->GetProperty(0);
        if (verboseLevel > 0) {
            G4cout << " YieldRatio = "<< YieldRatio << " and dump yield ratio table (yield ratio = fast/(fast+slow): " << G4endl;
            const_cast <G4MaterialPropertyVector*>(ptable)->DumpVector();
          }
    }


    //loop over fast/slow scintillations
    for (G4int scnt = 1; scnt <= nscnt; scnt++) {

        G4double ScintillationTime = 0.*ns;
        G4PhysicsOrderedFreeVector* ScintillationIntegral = NULL;
         
        if (scnt == 1) {//fast
            if (nscnt == 1) {
                if(Fast_Intensity){
                    ScintillationTime   = fastTimeConstant;
                    ScintillationIntegral =
                        (G4PhysicsOrderedFreeVector*)((*theFastIntegralTable)(materialIndex));
                }
                if(Slow_Intensity){
                    ScintillationTime   = slowTimeConstant;
                    ScintillationIntegral =
                        (G4PhysicsOrderedFreeVector*)((*theSlowIntegralTable)(materialIndex));
                }
            }
            else {
                if ( ExcitationRatio == 1.0 ) {
                  Num = G4int( 0.5 +  (min(YieldRatio,1.0) * NumTracks) );  // round off, not truncation
                }
                else {
                  Num = G4int( 0.5 +  (min(ExcitationRatio,1.0) * NumTracks));
                }
                if ( verboseLevel>1 ){
                  G4cout << "Generate Num " << Num << " optical photons with fast component using NumTracks " 
                         << NumTracks << " YieldRatio " << YieldRatio << " ExcitationRatio " << ExcitationRatio 
                         << " min(YieldRatio,1.)*NumTracks = " <<  min(YieldRatio,1.)*NumTracks 
                         << " min(ExcitationRatio,1.)*NumTracks = " <<  min(ExcitationRatio,1.)*NumTracks 
                         << G4endl;
                }
                ScintillationTime   = fastTimeConstant;
                ScintillationIntegral =
                    (G4PhysicsOrderedFreeVector*)((*theFastIntegralTable)(materialIndex));
            }
        }
        else {//slow
            Num = NumTracks - Num;
            ScintillationTime   =   slowTimeConstant;
            ScintillationIntegral =
                (G4PhysicsOrderedFreeVector*)((*theSlowIntegralTable)(materialIndex));
        }
        if (verboseLevel > 0) {
          G4cout << "generate " << Num << " optical photons with scintTime " << ScintillationTime 
                 << " slowTimeConstant " << slowTimeConstant << " fastTimeConstant " << fastTimeConstant << G4endl;
        }

        if (!ScintillationIntegral) continue;
        
        // Max Scintillation Integral
                
        for (G4int i = 0; i < Num; i++) { //Num is # of 2ndary tracks now
            // Determine photon energy

        if(scnt == 2) {
            ScintillationTime   =   slowTimeConstant;
            if((G4UniformRand() < slowerRatio) && (!flagReemission)) { 
              ScintillationTime = slowerTimeConstant;
            }
        }

            G4double sampledEnergy;
            if ( !flagReemission ) {
                // normal scintillation
                G4double CIIvalue = G4UniformRand()*
                    ScintillationIntegral->GetMaxValue();
                sampledEnergy=
                    ScintillationIntegral->GetEnergy(CIIvalue);

                if (verboseLevel>1) 
                    {
                        G4cout << "sampledEnergy = " << sampledEnergy << G4endl;
                        G4cout << "CIIvalue =        " << CIIvalue << G4endl;
                    }
            }
            else {
                // reemission, the sample method need modification
                G4double CIIvalue = G4UniformRand()*
                    ScintillationIntegral->GetMaxValue();
                if (CIIvalue == 0.0) {
                    // return unchanged particle and no secondaries  
                    aParticleChange.SetNumberOfSecondaries(0);
                    return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
                }
                sampledEnergy=
                    ScintillationIntegral->GetEnergy(CIIvalue);
                if (verboseLevel>1) {
                    G4cout << "oldEnergy = " <<aTrack.GetKineticEnergy() << G4endl;
                    G4cout << "reemittedSampledEnergy = " << sampledEnergy
                           << "\nreemittedCIIvalue =        " << CIIvalue << G4endl;
                }
            }

            // Generate random photon direction

            G4double cost = 1. - 2.*G4UniformRand();
            G4double sint = sqrt((1.-cost)*(1.+cost));

            G4double phi = twopi*G4UniformRand();
            G4double sinp = sin(phi);
            G4double cosp = 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 = sin(phi);
            cosp = cos(phi);

            photonPolarization = cosp * photonPolarization + sinp * perp;

            photonPolarization = photonPolarization.unit();

            // Generate a new photon:

            G4DynamicParticle* aScintillationPhoton =
                new G4DynamicParticle(G4OpticalPhoton::OpticalPhoton(), 
                                      photonMomentum);
            aScintillationPhoton->SetPolarization
                (photonPolarization.x(),
                 photonPolarization.y(),
                 photonPolarization.z());

            aScintillationPhoton->SetKineticEnergy(sampledEnergy);

            // Generate new G4Track object:

            G4double rand=0;
            G4ThreeVector aSecondaryPosition;
            G4double deltaTime;
            if (flagReemission) {
                deltaTime= pPostStepPoint->GetGlobalTime() - t0
                           -ScintillationTime * log( G4UniformRand() );
                aSecondaryPosition= pPostStepPoint->GetPosition();
                vertpos = aTrack.GetVertexPosition();
                vertenergy = aTrack.GetKineticEnergy();
                reem_d = 
                    sqrt( pow( aSecondaryPosition.x()-vertpos.x(), 2)
                          + pow( aSecondaryPosition.y()-vertpos.y(), 2)
                          + pow( aSecondaryPosition.z()-vertpos.z(), 2) );
            }
            else {
                if (aParticle->GetDefinition()->GetPDGCharge() != 0) 
                    {
                        rand = G4UniformRand();
                    }
                else
                    {
                        rand = 1.0;
                    }

                G4double delta = rand * aStep.GetStepLength();
                deltaTime = delta /
                    ((pPreStepPoint->GetVelocity()+
                      pPostStepPoint->GetVelocity())/2.);

                deltaTime = deltaTime - 
                    ScintillationTime * log( G4UniformRand() );

                aSecondaryPosition =
                    x0 + rand * aStep.GetDeltaPosition();
            }
            G4double aSecondaryTime = t0 + deltaTime;

            if ( verboseLevel>1 ){
              G4cout << "Generate " << i << "th scintillation photon at relative time(ns) " << deltaTime 
                     << " with ScintillationTime " << ScintillationTime << " flagReemission " << flagReemission << G4endl;
            }
            G4Track* aSecondaryTrack = 
                new G4Track(aScintillationPhoton,aSecondaryTime,aSecondaryPosition);

            G4CompositeTrackInfo* comp=new G4CompositeTrackInfo();
            DsPhotonTrackInfo* trackinf=new DsPhotonTrackInfo();
            if ( flagReemission ){
                if ( reemittedTI ) *trackinf = *reemittedTI;
                trackinf->SetReemitted();
            }
            else if ( fApplyPreQE ) {
                trackinf->SetMode(DsPhotonTrackInfo::kQEPreScale);
                trackinf->SetQE(fPreQE);
            }
            comp->SetPhotonTrackInfo(trackinf);
            aSecondaryTrack->SetUserInformation(comp);
                
            aSecondaryTrack->SetWeight( weight );
            aSecondaryTrack->SetTouchableHandle(aStep.GetPreStepPoint()->GetTouchableHandle());
            // aSecondaryTrack->SetTouchableHandle((G4VTouchable*)0);//this is wrong
                
            aSecondaryTrack->SetParentID(aTrack.GetTrackID());
                
            // add the secondary to the ParticleChange object
            aParticleChange.SetSecondaryWeightByProcess( true ); // recommended
            aParticleChange.AddSecondary(aSecondaryTrack);
                
            aSecondaryTrack->SetWeight( weight );
            if ( verboseLevel > 0 ) {
              G4cout << " aSecondaryTrack->SetWeight( " << weight<< " ) ; aSecondaryTrack->GetWeight() = " << aSecondaryTrack->GetWeight() << G4endl;}
            // The above line is necessary because AddSecondary() 
            // overrides our setting of the secondary track weight, 
            // in Geant4.3.1 & earlier. (and also later, at least 
            // until Geant4.7 (and beyond?)
            //  -- maybe not required if SetWeightByProcess(true) called,
            //  but we do both, just to be sure)
        }
    } // end loop over fast/slow scints

    if (verboseLevel > 0) {
        G4cout << "\n Exiting from G4Scintillation::DoIt -- NumberOfSecondaries = " 
               << aParticleChange.GetNumberOfSecondaries() << G4endl;
    }

    return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
}
G4VParticleChange * DsG4Scintillation::AtRestDoIt ( const G4Track &  aTrack,
const G4Step &  aStep 
)

Definition at line 158 of file DsG4Scintillation.cc.

{
    return DsG4Scintillation::PostStepDoIt(aTrack, aStep);
}
void DsG4Scintillation::SetTrackSecondariesFirst ( const G4bool  state) [inline]

Definition at line 311 of file DsG4Scintillation.h.

G4bool DsG4Scintillation::GetTrackSecondariesFirst ( ) const [inline]

Definition at line 317 of file DsG4Scintillation.h.

void DsG4Scintillation::SetScintillationYieldFactor ( const G4double  yieldfactor) [inline]

Definition at line 323 of file DsG4Scintillation.h.

{
        YieldFactor = yieldfactor;
}
G4double DsG4Scintillation::GetScintillationYieldFactor ( ) const [inline]

Definition at line 330 of file DsG4Scintillation.h.

{
        return YieldFactor;
}
void DsG4Scintillation::SetUseFastMu300nsTrick ( const G4bool  fastMu300nsTrick) [inline]

Definition at line 336 of file DsG4Scintillation.h.

{
        FastMu300nsTrick = fastMu300nsTrick;
}
G4bool DsG4Scintillation::GetUseFastMu300nsTrick ( ) const [inline]

Definition at line 342 of file DsG4Scintillation.h.

{
      return  FastMu300nsTrick;
}
void DsG4Scintillation::SetScintillationExcitationRatio ( const G4double  excitationratio) [inline]

Definition at line 351 of file DsG4Scintillation.h.

{
        ExcitationRatio = excitationratio;
}
G4double DsG4Scintillation::GetScintillationExcitationRatio ( ) const [inline]

Definition at line 357 of file DsG4Scintillation.h.

{
        return ExcitationRatio;
}
G4PhysicsTable * DsG4Scintillation::GetFastIntegralTable ( ) const [inline]

Definition at line 369 of file DsG4Scintillation.h.

G4PhysicsTable * DsG4Scintillation::GetSlowIntegralTable ( ) const [inline]

Definition at line 363 of file DsG4Scintillation.h.

G4PhysicsTable * DsG4Scintillation::GetReemissionIntegralTable ( ) const [inline]

Definition at line 375 of file DsG4Scintillation.h.

void DsG4Scintillation::DumpPhysicsTable ( ) const [inline]

Definition at line 381 of file DsG4Scintillation.h.

{
        if (theFastIntegralTable) {
           G4int PhysicsTableSize = theFastIntegralTable->entries();
           G4PhysicsOrderedFreeVector *v;

           for (G4int i = 0 ; i < PhysicsTableSize ; i++ )
           {
                v = (G4PhysicsOrderedFreeVector*)(*theFastIntegralTable)[i];
                v->DumpValues();
           }
         }

        if (theSlowIntegralTable) {
           G4int PhysicsTableSize = theSlowIntegralTable->entries();
           G4PhysicsOrderedFreeVector *v;

           for (G4int i = 0 ; i < PhysicsTableSize ; i++ )
           {
                v = (G4PhysicsOrderedFreeVector*)(*theSlowIntegralTable)[i];
                v->DumpValues();
           }
         }

        if (theReemissionIntegralTable) {
           G4int PhysicsTableSize = theReemissionIntegralTable->entries();
           G4PhysicsOrderedFreeVector *v;

           for (G4int i = 0 ; i < PhysicsTableSize ; i++ )
           {
                v = (G4PhysicsOrderedFreeVector*)(*theReemissionIntegralTable)[i];
                v->DumpValues();
           }
         }
}
G4double DsG4Scintillation::GetPhotonWeight ( ) const [inline]

Definition at line 202 of file DsG4Scintillation.h.

{ return fPhotonWeight; }
void DsG4Scintillation::SetPhotonWeight ( G4double  weight) [inline]

Definition at line 203 of file DsG4Scintillation.h.

{ fPhotonWeight = weight; }
void DsG4Scintillation::SetDoReemission ( bool  tf = true) [inline]

Definition at line 204 of file DsG4Scintillation.h.

{ doReemission = tf; }
bool DsG4Scintillation::GetDoReemission ( ) [inline]

Definition at line 205 of file DsG4Scintillation.h.

{ return doReemission; }
void DsG4Scintillation::SetDoBothProcess ( bool  tf = true) [inline]

Definition at line 206 of file DsG4Scintillation.h.

{ doBothProcess = tf; }
bool DsG4Scintillation::GetDoBothProcess ( ) [inline]

Definition at line 207 of file DsG4Scintillation.h.

{ return doBothProcess; }
G4bool DsG4Scintillation::GetApplyPreQE ( ) const [inline]

Definition at line 209 of file DsG4Scintillation.h.

{ return fApplyPreQE; }
void DsG4Scintillation::SetApplyPreQE ( G4bool  a) [inline]

Definition at line 210 of file DsG4Scintillation.h.

{ fApplyPreQE = a; }
G4double DsG4Scintillation::GetPreQE ( ) const [inline]

Definition at line 212 of file DsG4Scintillation.h.

{ return fPreQE; }
void DsG4Scintillation::SetPreQE ( G4double  a) [inline]

Definition at line 213 of file DsG4Scintillation.h.

{ fPreQE = a; }
void DsG4Scintillation::SetBirksConstant1 ( double  c1) [inline]

Definition at line 215 of file DsG4Scintillation.h.

double DsG4Scintillation::GetBirksConstant1 ( ) [inline]

Definition at line 216 of file DsG4Scintillation.h.

{return birksConstant1;}
void DsG4Scintillation::SetBirksConstant2 ( double  c2) [inline]

Definition at line 217 of file DsG4Scintillation.h.

double DsG4Scintillation::GetBirksConstant2 ( ) [inline]

Definition at line 218 of file DsG4Scintillation.h.

{return birksConstant2;}
void DsG4Scintillation::SetGammaSlowerTimeConstant ( double  st) [inline]

Definition at line 220 of file DsG4Scintillation.h.

{ gammaSlowerTime = st;}
double DsG4Scintillation::GetGammaSlowerTimeConstant ( ) [inline]

Definition at line 221 of file DsG4Scintillation.h.

{return gammaSlowerTime;}
void DsG4Scintillation::SetGammaSlowerRatio ( double  sr) [inline]

Definition at line 223 of file DsG4Scintillation.h.

double DsG4Scintillation::GetGammaSlowerRatio ( ) [inline]

Definition at line 224 of file DsG4Scintillation.h.

{return gammaSlowerRatio;}
void DsG4Scintillation::SetNeutronSlowerTimeConstant ( double  st) [inline]

Definition at line 226 of file DsG4Scintillation.h.

double DsG4Scintillation::GetNeutronSlowerTimeConstant ( ) [inline]

Definition at line 227 of file DsG4Scintillation.h.

void DsG4Scintillation::SetNeutronSlowerRatio ( double  sr) [inline]

Definition at line 229 of file DsG4Scintillation.h.

double DsG4Scintillation::GetNeutronSlowerRatio ( ) [inline]

Definition at line 230 of file DsG4Scintillation.h.

void DsG4Scintillation::SetAlphaSlowerTimeConstant ( double  st) [inline]

Definition at line 232 of file DsG4Scintillation.h.

{ alphaSlowerTime = st;}
double DsG4Scintillation::GetAlphaSlowerTimeConstant ( ) [inline]

Definition at line 233 of file DsG4Scintillation.h.

{return alphaSlowerTime;}
void DsG4Scintillation::SetAlphaSlowerRatio ( double  sr) [inline]

Definition at line 235 of file DsG4Scintillation.h.

double DsG4Scintillation::GetAlphaSlowerRatio ( ) [inline]

Definition at line 236 of file DsG4Scintillation.h.

{return alphaSlowerRatio;}
void DsG4Scintillation::SetNoOp ( bool  tf = true) [inline]

Definition at line 243 of file DsG4Scintillation.h.

{ m_noop = tf; }
void DsG4Scintillation::BuildThePhysicsTable ( ) [private]

Definition at line 735 of file DsG4Scintillation.cc.

{
    if (theFastIntegralTable && theSlowIntegralTable && theReemissionIntegralTable) return;

    const G4MaterialTable* theMaterialTable = 
        G4Material::GetMaterialTable();
    G4int numOfMaterials = G4Material::GetNumberOfMaterials();

    // create new physics table
    if (verboseLevel > 0) {
      G4cout << " theFastIntegralTable " << theFastIntegralTable 
             << " theSlowIntegralTable " << theSlowIntegralTable 
             << " theReemissionIntegralTable " << theReemissionIntegralTable << G4endl;
    }
    if(!theFastIntegralTable)theFastIntegralTable = new G4PhysicsTable(numOfMaterials);
    if(!theSlowIntegralTable)theSlowIntegralTable = new G4PhysicsTable(numOfMaterials);
    if(!theReemissionIntegralTable)theReemissionIntegralTable
                                       = new G4PhysicsTable(numOfMaterials);
    if (verboseLevel > 0) {
      G4cout << " building the physics tables for the scintillation process " << G4endl;
    }
    // loop for materials

    for (G4int i=0 ; i < numOfMaterials; i++) {
        G4PhysicsOrderedFreeVector* aPhysicsOrderedFreeVector =
            new G4PhysicsOrderedFreeVector();
        G4PhysicsOrderedFreeVector* bPhysicsOrderedFreeVector =
            new G4PhysicsOrderedFreeVector();
        G4PhysicsOrderedFreeVector* cPhysicsOrderedFreeVector =
            new G4PhysicsOrderedFreeVector();

        // 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) {
              if (verboseLevel > 0) {
                G4cout << " Building the material properties table for FASTCOMPONENT" << G4endl;
              }
                // 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;
                    }

                }
            }

            G4MaterialPropertyVector* theSlowLightVector =
                aMaterialPropertiesTable->GetProperty("SLOWCOMPONENT");

            if (theSlowLightVector) {
              if (verboseLevel > 0) {
                G4cout << " Building the material properties table for SLOWCOMPONENT" << G4endl;
              }

                // Retrieve the first intensity point in vector
                // of (photon energy, intensity) pairs

                theSlowLightVector->ResetIterator();
                ++(*theSlowLightVector);  // advance to 1st entry

                G4double currentIN = theSlowLightVector->
                    GetProperty();

                if (currentIN >= 0.0) {

                    // Create first (photon energy, Scintillation
                    // Integral pair

                    G4double currentPM = theSlowLightVector->
                        GetPhotonEnergy();

                    G4double currentCII = 0.0;

                    bPhysicsOrderedFreeVector->
                        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(++(*theSlowLightVector)) {
                        currentPM = theSlowLightVector->
                            GetPhotonEnergy();

                        currentIN=theSlowLightVector->
                            GetProperty();

                        currentCII = 0.5 * (prevIN + currentIN);

                        currentCII = prevCII +
                            (currentPM - prevPM) * currentCII;

                        bPhysicsOrderedFreeVector->
                            InsertValues(currentPM, currentCII);

                        prevPM  = currentPM;
                        prevCII = currentCII;
                        prevIN  = currentIN;
                    }

                }
            }

            G4MaterialPropertyVector* theReemissionVector =
                aMaterialPropertiesTable->GetProperty("REEMISSIONPROB");

            if (theReemissionVector) {
              if (verboseLevel > 0) {
                G4cout << " Building the material properties table for REEMISSIONPROB" << G4endl;
              }

                // Retrieve the first intensity point in vector
                // of (photon energy, intensity) pairs

                theReemissionVector->ResetIterator();
                ++(*theReemissionVector);  // advance to 1st entry

                G4double currentIN = theReemissionVector->
                    GetProperty();

                if (currentIN >= 0.0) {

                    // Create first (photon energy, Scintillation
                    // Integral pair

                    G4double currentPM = theReemissionVector->
                        GetPhotonEnergy();

                    G4double currentCII = 0.0;

                    cPhysicsOrderedFreeVector->
                        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(++(*theReemissionVector)) {
                        currentPM = theReemissionVector->
                            GetPhotonEnergy();

                        currentIN=theReemissionVector->
                            GetProperty();

                        currentCII = 0.5 * (prevIN + currentIN);

                        currentCII = prevCII +
                            (currentPM - prevPM) * currentCII;

                        cPhysicsOrderedFreeVector->
                            InsertValues(currentPM, currentCII);

                        prevPM  = currentPM;
                        prevCII = currentCII;
                        prevIN  = currentIN;
                    }

                }
            }

        }

        // 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);
        theReemissionIntegralTable->insertAt(i,cPhysicsOrderedFreeVector);
    }
}

Member Data Documentation

G4PhysicsTable* DsG4Scintillation::theSlowIntegralTable [protected]

Definition at line 256 of file DsG4Scintillation.h.

G4PhysicsTable* DsG4Scintillation::theFastIntegralTable [protected]

Definition at line 257 of file DsG4Scintillation.h.

G4PhysicsTable* DsG4Scintillation::theReemissionIntegralTable [protected]

Definition at line 258 of file DsG4Scintillation.h.

G4bool DsG4Scintillation::doReemission [protected]

Definition at line 261 of file DsG4Scintillation.h.

Definition at line 263 of file DsG4Scintillation.h.

Definition at line 266 of file DsG4Scintillation.h.

Definition at line 267 of file DsG4Scintillation.h.

Definition at line 269 of file DsG4Scintillation.h.

double DsG4Scintillation::slowerRatio [protected]

Definition at line 270 of file DsG4Scintillation.h.

Definition at line 272 of file DsG4Scintillation.h.

Definition at line 273 of file DsG4Scintillation.h.

Definition at line 274 of file DsG4Scintillation.h.

Definition at line 275 of file DsG4Scintillation.h.

Definition at line 276 of file DsG4Scintillation.h.

Definition at line 277 of file DsG4Scintillation.h.

Definition at line 282 of file DsG4Scintillation.h.

G4double DsG4Scintillation::YieldFactor [private]

Definition at line 284 of file DsG4Scintillation.h.

Definition at line 285 of file DsG4Scintillation.h.

Definition at line 286 of file DsG4Scintillation.h.

Definition at line 289 of file DsG4Scintillation.h.

Definition at line 290 of file DsG4Scintillation.h.

G4double DsG4Scintillation::fPreQE [private]

Definition at line 291 of file DsG4Scintillation.h.

bool DsG4Scintillation::m_noop [private]

Definition at line 292 of file DsG4Scintillation.h.


The documentation for this class was generated from the following files:
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:17:58 for DetSim by doxygen 1.7.4