EDepSim::DokeBirks Class Reference

#include <EDepSimDokeBirks.hh>

Inheritance diagram for EDepSim::DokeBirks:

Public Member Functions
	DokeBirks (const G4String &processName="Doke-Birks", G4ProcessType type=fElectromagnetic)
	~DokeBirks ()
G4bool	IsApplicable (const G4ParticleDefinition &aParticleType)
	Determine which particles this process should be applied too. More...
G4double	GetMeanFreePath (const G4Track &aTrack, G4double, G4ForceCondition *)
G4double	GetMeanLifeTime (const G4Track &aTrack, G4ForceCondition *)
G4VParticleChange *	PostStepDoIt (const G4Track &aTrack, const G4Step &aStep)
	Apply the scintilation process for an in-flight particle. More...
G4VParticleChange *	AtRestDoIt (const G4Track &aTrack, const G4Step &aStep)

Private Member Functions
G4double	CalculateElectronLET (G4double E)

Private Attributes
G4EmSaturation *	fEmSaturation

Detailed Description

Definition at line 31 of file EDepSimDokeBirks.hh.

Constructor & Destructor Documentation

EDepSim::DokeBirks::DokeBirks	(	const G4String &	processName = "Doke-Birks",
		G4ProcessType	type = fElectromagnetic
	)

Definition at line 18 of file EDepSimDokeBirks.cc.

    : G4VRestDiscreteProcess(processName, type),
      fEmSaturation(nullptr)
{
    SetProcessSubType(fScintillation);
    
    if (verboseLevel>0) {
        G4cout << GetProcessName() << " is created " << G4endl;
    }

}

EDepSim::DokeBirks::~DokeBirks

(

)

Definition at line 31 of file EDepSimDokeBirks.cc.

{} //destructor needed to avoid linker error

Member Function Documentation

G4VParticleChange * EDepSim::DokeBirks::AtRestDoIt	(	const G4Track &	aTrack,
		const G4Step &	aStep
	)

Apply the scintilation process to a stopped particle (this just calles PostStepDoIt().

Definition at line 52 of file EDepSimDokeBirks.cc.

{
    return EDepSim::DokeBirks::PostStepDoIt(aTrack, aStep);
}

G4double EDepSim::DokeBirks::CalculateElectronLET ( G4double  E )

private

Get the expected electron electron dEdX as a function of energy. This is normalized to denstity. The energy must be in MeV.

Definition at line 233 of file EDepSimDokeBirks.cc.

                                                            {
    G4double LET;
    if ( E >= 1 ) LET = 116.70-162.97*log10(E)+99.361*pow(log10(E),2)-
                      33.405*pow(log10(E),3)+6.5069*pow(log10(E),4)-
                      0.69334*pow(log10(E),5)+.031563*pow(log10(E),6);
    else if ( E>0 && E<1 ) LET = 100;
    else LET = 0;
    return LET;
}

G4double EDepSim::DokeBirks::GetMeanFreePath	(	const G4Track &	aTrack,
		G4double	,
		G4ForceCondition *	condition
	)

Don't limit the step since this process doesn't actually change the energy deposit (it returns infinity). The process does set the 'StronglyForced' condition so the DoIt is invoked at every step.

Definition at line 206 of file EDepSimDokeBirks.cc.

{
    *condition = StronglyForced;
    // Force the EDepSim::DokeBirks physics process to always happen, even if it's
    // not the dominant process.  In effect it is a meta-process on top of any
    // and all other energy depositions which may occur, just like the
    // original G4Scintillation (disregard DBL_MAX, this function makes the
    // mean free path zero really, not infinite)

    return DBL_MAX;
}

G4double EDepSim::DokeBirks::GetMeanLifeTime	(	const G4Track &	aTrack,
		G4ForceCondition *	condition
	)

Don't limit the step since this process doesn't actually change the energy deposit (it returns infinity). The process does set the 'StronglyForced' condition so the DoIt is invoked at every step.

Definition at line 222 of file EDepSimDokeBirks.cc.

{
    *condition = Forced;
    // This function and this condition has the same effect as the above
    return DBL_MAX;
}

G4bool EDepSim::DokeBirks::IsApplicable

(

const G4ParticleDefinition &

aParticleType

)

Determine which particles this process should be applied too.

Definition at line 35 of file EDepSimDokeBirks.cc.

{
    // This process is applicable to some neutral particles (e.g. gammas)
    // since G4 will fake creating the low energy electrons.
    if (aParticleType.IsShortLived()) return false;
    if (aParticleType.GetParticleName() == "opticalphoton") return false;
    if (aParticleType.GetParticleName() == "thermalelectron") return false;
    // Not applicable to the neutrinos.
    if (std::abs(aParticleType.GetPDGEncoding())==12) return false;
    if (std::abs(aParticleType.GetPDGEncoding())==14) return false;
    if (std::abs(aParticleType.GetPDGEncoding())==16) return false;
    return true;
}

G4VParticleChange * EDepSim::DokeBirks::PostStepDoIt	(	const G4Track &	aTrack,
		const G4Step &	aStep
	)

Apply the scintilation process for an in-flight particle.

Definition at line 60 of file EDepSimDokeBirks.cc.

{
    aParticleChange.Initialize(aTrack);
        
    const G4Material* aMaterial = aTrack.GetMaterial();

    bool inLiquidArgon = true;
    // Check that we are in liquid argon.
    if (aMaterial->GetState() != kStateLiquid) {
        // Do Nothing
        if (aMaterial->GetState() == kStateUndefined) {
            EDepSimError("Undefined material state for "
                         << aMaterial->GetName());
        }
        inLiquidArgon = false;
    }
    if (aMaterial->GetNumberOfElements() > 1) {
        // Do Nothing
        inLiquidArgon = false;
    }
    if (std::abs(aMaterial->GetElement(0)->GetZ()-18.0) > 0.5) {
        // Do Nothing
        inLiquidArgon = false;
    }
    
    const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
    const G4ParticleDefinition *pDef = aParticle->GetDefinition();
    G4String particleName = pDef->GetParticleName();
    // G4StepPoint* pPreStepPoint  = aStep.GetPreStepPoint();
    G4double totalEnergyDeposit = aStep.GetTotalEnergyDeposit();
    
    if (totalEnergyDeposit <= 0.0) {
        return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);        
    }

    if (aStep.GetStepLength() <= 0.0) {
        return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);        
    }
    
    if (!inLiquidArgon) {
        if (!fEmSaturation) fEmSaturation = new G4EmSaturation(0);
        G4double nonIonizingEnergy
            = fEmSaturation->VisibleEnergyDepositionAtAStep(&aStep);
        aParticleChange.ProposeNonIonizingEnergyDeposit(nonIonizingEnergy);
        return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
    }

    // Figure out the electric field for the volume.
    G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
    const G4VPhysicalVolume* aVolume = pPostStepPoint->GetPhysicalVolume();
    const G4LogicalVolume* aLogVolume = aVolume->GetLogicalVolume();
    const G4FieldManager* aFieldManager = aLogVolume->GetFieldManager();
    if (!aFieldManager) {
        return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
    }
    if (!aFieldManager->DoesFieldExist()) {
        return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
    }
    const G4Field* aField = aFieldManager->GetDetectorField();
    if (!aField) {
        EDepSimError("LAr Volume "
                     << aLogVolume->GetName()
                     << " should have field!");
        return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);        
    }

    G4ThreeVector thePosition = pPostStepPoint->GetPosition();
    double theTime = pPostStepPoint->GetGlobalTime();
    double point[4] = {
        thePosition.x(), thePosition.y(), thePosition.z(), theTime};
    double bField[6];
    aField->GetFieldValue(point,bField);
    
    double electricField = bField[3]*bField[3];
    electricField += bField[4]*bField[4];
    electricField += bField[5]*bField[5];
    if (electricField > 0) {
        electricField = std::sqrt(electricField);
    }
    // The code below is pulled from G4S1Light and is simplified to be
    // Doke-Birks only, in LAr only, and for an electric field only.  This is
    // for ARGON only.  The Doke-Birks constants are in kilovolt/cm
    G4double dokeBirks[3];

    dokeBirks[0] = 0.07*pow((electricField),-0.85);
    dokeBirks[2] = 0.00;
    dokeBirks[1] = dokeBirks[0]/(1-dokeBirks[2]); //B=A/(1-C) (see paper)

    G4double dE = totalEnergyDeposit/MeV;
    G4double dx = aStep.GetStepLength()/cm;
    G4double density = aMaterial->GetDensity()/(g/cm3);
    G4double LET = 0.0; //lin. energy transfer (prop. to dE/dx)
    if (dx != 0.0) LET = (dE/dx)*(1/density);

    // There was some dEdX for a photon.  This means that the very low energy
    // electrons were not simulated.
    if (particleName == "gamma") {
        LET = CalculateElectronLET( 1000*dE );
        dx = dE/density/LET;
    }
    else if (std::abs(pDef->GetPDGCharge()) < 0.1) {
        // Handle the case where a neutral particle made an "electron".  This
        // is when there is a tiny electron scatter that is below the G4
        // threshold.  The assumption is the same as NEST.  The product
        // particle is an electron.
        double tempLET = CalculateElectronLET( 1000*dE );
        double electronStep = dE/density/tempLET;
        if (dx > electronStep) {
            LET = tempLET;
            dx = electronStep;
        }
    }

    // Special case for electrons that G4 may stop when they get to short and
    // truncate the energy deposition.
    do {
        if (particleName != "e+" || particleName != "e-") break;
        if (aTrack.GetCurrentStepNumber() != 1) break;
        double ratio = CalculateElectronLET( 1000.0*dE ) / LET;
        if (ratio > 0.7) break;
        LET *= ratio;
        dx /= ratio;
    } while (false);
    
    G4double recombProb = (dokeBirks[0]*LET)/(1+dokeBirks[1]*LET)+dokeBirks[2];

    G4double ScintillationYield = 1 / (19.5*eV);
    G4double NumQuanta = ScintillationYield*totalEnergyDeposit;
    G4double ExcitationRatio = 0.21; //Aprile et. al book
    G4double NumExcitons = NumQuanta*ExcitationRatio/(1+ExcitationRatio);
    G4double NumIons = NumQuanta - NumExcitons;
    G4double NumPhotons = NumExcitons + NumIons*recombProb;

    // The number of photons and electrons for this step has been calculated
    // using the Doke/Birks law version of recombination.  For DETSIM we are
    // only looking at argon, and according the Matt both Doke-Birks and
    // Thomas-Imel fit the data.  It's also because we are not interested in
    // extremely low energies, so Thomas-Imel makes a negligible correction.
    G4double nonIonizingEnergy = totalEnergyDeposit;
    nonIonizingEnergy *= NumPhotons/NumQuanta; 
    aParticleChange.ProposeNonIonizingEnergyDeposit(nonIonizingEnergy);
    return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
}

Member Data Documentation

G4EmSaturation* EDepSim::DokeBirks::fEmSaturation

private

Definition at line 71 of file EDepSimDokeBirks.hh.

---

The documentation for this class was generated from the following files:

• edep-sim/src/EDepSimDokeBirks.hh
• edep-sim/src/EDepSimDokeBirks.cc