G4DNAPTBExcitationModel.cc

Go to the documentation of this file.

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
// Authors: S. Meylan and C. Villagrasa (IRSN, France)
// Models come from
// M. Bug et al, Rad. Phys and Chem. 130, 459-479 (2017)
//
 
#include "G4DNAPTBExcitationModel.hh"
#include "G4SystemOfUnits.hh"
#include "G4DNAChemistryManager.hh"
#include "G4DNAMolecularMaterial.hh"
 
G4DNAPTBExcitationModel::G4DNAPTBExcitationModel(const G4String& applyToMaterial, const G4ParticleDefinition*,
                                                   const G4String& nam)
    : G4VDNAModel(nam, applyToMaterial)
{
    verboseLevel= 0;
    // Verbosity scale:
    // 0 = nothing
    // 1 = warning for energy non-conservation
    // 2 = details of energy budget
    // 3 = calculation of cross sections, file openings, sampling of atoms
    // 4 = entering in methods
 
    // initialisation of mean energy loss for each material
    tableMeanEnergyPTB["THF"] = 8.01*eV;
    tableMeanEnergyPTB["PY"] = 7.61*eV;
    tableMeanEnergyPTB["PU"] = 7.61*eV;
    tableMeanEnergyPTB["TMP"] = 8.01*eV;
 
    if( verboseLevel>0 )
    {
        G4cout << "PTB excitation model is constructed " << G4endl;
    }
}
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
G4DNAPTBExcitationModel::~G4DNAPTBExcitationModel()
{ 
 
}
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
void G4DNAPTBExcitationModel::Initialise(const G4ParticleDefinition* particle,
                                         const G4DataVector& /*cuts*/, G4ParticleChangeForGamma*)
{
    if (verboseLevel > 3)
        G4cout << "Calling G4DNAPTBExcitationModel::Initialise()" << G4endl;
 
    G4double scaleFactor = 1e-16*cm*cm;
    G4double scaleFactorBorn = (1.e-22 / 3.343) * m*m;
 
    G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
 
    //*******************************************************
    // Cross section data
    //*******************************************************
 
    if(particle == electronDef)
    {
        G4String particleName = particle->GetParticleName();
 
        AddCrossSectionData("THF",
                            particleName,
                            "dna/sigma_excitation_e-_PTB_THF",
                            scaleFactor);
        SetLowELimit("THF", particleName, 9.*eV);
        SetHighELimit("THF", particleName, 1.*keV);
 
        AddCrossSectionData("PY",
                            particleName,
                            "dna/sigma_excitation_e-_PTB_PY",
                            scaleFactor);
        SetLowELimit("PY", particleName, 9.*eV);
        SetHighELimit("PY", particleName, 1.*keV);
 
        AddCrossSectionData("PU",
                            particleName,
                            "dna/sigma_excitation_e-_PTB_PU",
                            scaleFactor);
        SetLowELimit("PU", particleName, 9.*eV);
        SetHighELimit("PU", particleName, 1.*keV);
 
        AddCrossSectionData("TMP",
                            particleName,
                            "dna/sigma_excitation_e-_PTB_TMP",
                            scaleFactor);
        SetLowELimit("TMP", particleName, 9.*eV);
        SetHighELimit("TMP", particleName, 1.*keV);
 
        AddCrossSectionData("G4_WATER",
                            particleName,
                            "dna/sigma_excitation_e_born",
                            scaleFactorBorn);
        SetLowELimit("G4_WATER", particleName, 9.*eV);
        SetHighELimit("G4_WATER", particleName, 1.*keV);
 
        // DNA materials
        //
        AddCrossSectionData("backbone_THF",
                            particleName,
                            "dna/sigma_excitation_e-_PTB_THF",
                            scaleFactor*33./30);
        SetLowELimit("backbone_THF", particleName, 9.*eV);
        SetHighELimit("backbone_THF", particleName, 1.*keV);
 
        AddCrossSectionData("cytosine_PY",
                            particleName,
                            "dna/sigma_excitation_e-_PTB_PY",
                            scaleFactor*42./30);
        SetLowELimit("cytosine_PY", particleName, 9.*eV);
        SetHighELimit("cytosine_PY", particleName, 1.*keV);
 
        AddCrossSectionData("thymine_PY",
                            particleName,
                            "dna/sigma_excitation_e-_PTB_PY",
                            scaleFactor*48./30);
        SetLowELimit("thymine_PY", particleName, 9.*eV);
        SetHighELimit("thymine_PY", particleName, 1.*keV);
 
        AddCrossSectionData("adenine_PU",
                            particleName,
                            "dna/sigma_excitation_e-_PTB_PU",
                            scaleFactor*50./44);
        SetLowELimit("adenine_PU", particleName, 9.*eV);
        SetHighELimit("adenine_PU", particleName, 1.*keV);
 
        AddCrossSectionData("guanine_PU",
                            particleName,
                            "dna/sigma_excitation_e-_PTB_PU",
                            scaleFactor*56./44);
        SetLowELimit("guanine_PU", particleName, 9.*eV);
        SetHighELimit("guanine_PU", particleName, 1.*keV);
 
        AddCrossSectionData("backbone_TMP",
                            particleName,
                            "dna/sigma_excitation_e-_PTB_TMP",
                            scaleFactor*33./50);
        SetLowELimit("backbone_TMP", particleName, 9.*eV);
        SetHighELimit("backbone_TMP", particleName, 1.*keV);
 
        // MPietrzak, adding paths for N2
        AddCrossSectionData("N2",
            particleName,
            "dna/sigma_excitation_e-_PTB_N2",
            scaleFactor);
        SetLowELimit("N2", particleName, 13.*eV);
        SetHighELimit("N2", particleName, 1.02*MeV);
        // MPietrzak
    }
 
    //*******************************************************
    // Load data
    //*******************************************************
 
    LoadCrossSectionData(particle->GetParticleName() );
 
    //*******************************************************
    // Verbose
    //*******************************************************
 
    if( verboseLevel>0 )
    {
        G4cout << "PTB excitation model is initialized " << G4endl;
    }
}
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
G4double G4DNAPTBExcitationModel::CrossSectionPerVolume(const G4Material* /*material*/,
                                                        const G4String& materialName,
                                                        const G4ParticleDefinition* particleDefinition,
                                                        G4double ekin,
                                                        G4double /*emin*/,
                                                        G4double /*emax*/)
{
    if (verboseLevel > 3)
        G4cout << "Calling CrossSectionPerVolume() of G4DNAPTBExcitationModel" << G4endl;
 
    // Get the name of the current particle
    G4String particleName = particleDefinition->GetParticleName();
 
    // initialise variables
    G4double lowLim = 0;
    G4double highLim = 0;
    G4double sigma=0;
 
    // Get the low energy limit for the current particle
    lowLim = GetLowELimit(materialName, particleName);
 
    // Get the high energy limit for the current particle
    highLim = GetHighELimit(materialName, particleName);
 
    // Check that we are in the correct energy range
    if (ekin >= lowLim && ekin < highLim)
    {
        // Get the map with all the data tables
        TableMapData* tableData = GetTableData();
 
        // Retrieve the cross section value
        sigma = (*tableData)[materialName][particleName]->FindValue(ekin);
 
        if (verboseLevel > 2)
        {
            G4cout << "__________________________________" << G4endl;
            G4cout << "°°° G4DNAPTBExcitationModel - XS INFO START" << G4endl;
            G4cout << "°°° Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
            G4cout << "°°° Cross section per "<< materialName <<" molecule (cm^2)=" << sigma/cm/cm << G4endl;
            G4cout << "°°° G4DNAPTBExcitationModel - XS INFO END" << G4endl;
        }
 
    }
 
    // Return the cross section value
    return sigma;
}
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
void G4DNAPTBExcitationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
                                                const G4MaterialCutsCouple* /*couple*/,
                                                const G4String& materialName,
                                                const G4DynamicParticle* aDynamicParticle,
                                                G4ParticleChangeForGamma* particleChangeForGamma,
                                                G4double /*tmin*/,
                                                G4double /*tmax*/)
{
    if (verboseLevel > 3)
        G4cout << "Calling SampleSecondaries() of G4DNAPTBExcitationModel" << G4endl;
 
    // Get the incident particle kinetic energy
    G4double k = aDynamicParticle->GetKineticEnergy();
    //Get the particle name
    const G4String& particleName = aDynamicParticle->GetDefinition()->GetParticleName();
    // Get the energy limits
    G4double lowLim = GetLowELimit(materialName, particleName);
    G4double highLim = GetHighELimit(materialName, particleName);
 
    // Check if we are in the correct energy range
    if (k >= lowLim && k < highLim)
    {
        if(materialName=="N2")
        {
 
            // Retrieve the excitation energy for the current material
            G4int level = RandomSelectShell(k,particleName,materialName);
            G4double excitationEnergy = ptbExcitationStructure.ExcitationEnergy(level, materialName);
 
            // Calculate the new energy of the particle
            G4double newEnergy = k - excitationEnergy;
 
            // Check that the new energy is above zero before applying it the particle.
            // Otherwise, do nothing.
            if (newEnergy > 0)
            {
                particleChangeForGamma->ProposeMomentumDirection(aDynamicParticle->GetMomentumDirection());
                particleChangeForGamma->SetProposedKineticEnergy(newEnergy);
                particleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
                G4double ioniThres = ptbIonisationStructure.IonisationEnergy(0,materialName);
                // if excitation energy greater than ionisation threshold, then autoionisaiton
                if((excitationEnergy>ioniThres)&&(G4UniformRand()<0.5))
                {
                    particleChangeForGamma->ProposeLocalEnergyDeposit(ioniThres);
                    // energy of ejected electron
                    G4double secondaryKinetic = excitationEnergy - ioniThres;
                    // random direction
                    G4double cosTheta = 2*G4UniformRand() - 1., phi = CLHEP::twopi*G4UniformRand();
                    G4double sinTheta = std::sqrt(1. - cosTheta*cosTheta);
                    G4double ux = sinTheta*std::cos(phi),
                    uy = sinTheta*std::sin(phi),
                    uz = cosTheta;
                    G4ThreeVector deltaDirection(ux,uy,uz);
                    // Create the new particle with its characteristics
                    G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic) ;
                    fvect->push_back(dp);
                }
            } else {
                G4ExceptionDescription description;
                description<<"Kinetic energy <= 0 at "<<materialName<<" material !!!";
                G4Exception("G4DNAPTBExcitationModel::SampleSecondaries","",FatalException,description);
            }     
        } else if(materialName!="G4_WATER"){
            // Retrieve the excitation energy for the current material
            G4double excitationEnergy = tableMeanEnergyPTB[materialName];
            // Calculate the new energy of the particle
            G4double newEnergy = k - excitationEnergy;
            // Check that the new energy is above zero before applying it the particle.
            // Otherwise, do nothing.
            if (newEnergy > 0){
                particleChangeForGamma->ProposeMomentumDirection(aDynamicParticle->GetMomentumDirection());
                particleChangeForGamma->SetProposedKineticEnergy(newEnergy);
                particleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
            } else {
                G4ExceptionDescription description;
                description<<"Kinetic energy <= 0 at "<<materialName<<" material !!!";
                G4Exception("G4DNAPTBExcitationModel::SampleSecondaries","",FatalException,description);
            }
        } else {
            G4int level = RandomSelectShell(k,particleName, materialName);
            G4double excitationEnergy = waterStructure.ExcitationEnergy(level);
            G4double newEnergy = k - excitationEnergy;
 
            if (newEnergy > 0){
                particleChangeForGamma->ProposeMomentumDirection(aDynamicParticle->GetMomentumDirection());
                particleChangeForGamma->SetProposedKineticEnergy(newEnergy);
                particleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
                const G4Track * theIncomingTrack = particleChangeForGamma->GetCurrentTrack();
                G4DNAChemistryManager::Instance()->CreateWaterMolecule(eExcitedMolecule,
                                                                level,
                                                                theIncomingTrack);
            } else {
                G4ExceptionDescription description;
                description<<"Kinetic energy <= 0 at "<<materialName<<" material !!!";
                G4Exception("G4DNAPTBExcitationModel::SampleSecondaries","",FatalException,description);
            }
        }
    }
    
}