Ignore:
Timestamp:
Apr 6, 2009, 12:21:12 PM (15 years ago)
Author:
garnier
Message:

update processes

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/source/processes/electromagnetic/standard/src/G4ionIonisation.cc

    r819 r961  
    2424// ********************************************************************
    2525//
    26 // $Id: G4ionIonisation.cc,v 1.45.2.2 2008/04/25 00:34:55 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: G4ionIonisation.cc,v 1.66 2009/02/20 12:06:37 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
     
    5555// 16-05-07 Add data for light ion stopping only for GenericIon (V.Ivantchenko)
    5656// 07-11-07 Fill non-ionizing energy loss (V.Ivantchenko)
     57// 12-09-08 Removed InitialiseMassCharge and CorrectionsAlongStep (VI)
    5758//
    5859//
     
    6566#include "G4Electron.hh"
    6667#include "G4Proton.hh"
     68//#include "G4Alpha.hh"
    6769#include "G4GenericIon.hh"
    6870#include "G4BraggModel.hh"
    6971#include "G4BraggIonModel.hh"
    7072#include "G4BetheBlochModel.hh"
    71 #include "G4IonFluctuations.hh"
    7273#include "G4UnitsTable.hh"
    7374#include "G4LossTableManager.hh"
    7475#include "G4WaterStopping.hh"
     76#include "G4EmCorrections.hh"
     77#include "G4IonFluctuations.hh"
    7578
    7679//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     
    8083G4ionIonisation::G4ionIonisation(const G4String& name)
    8184  : G4VEnergyLossProcess(name),
     85    corr(0),
    8286    theParticle(0),
    83     theBaseParticle(0),
    8487    isInitialised(false),
    8588    stopDataActive(true),
     
    8992  SetStepFunction(0.1, 0.1*mm);
    9093  SetIntegral(true);
    91   SetVerboseLevel(1);
     94  SetProcessSubType(fIonisation);
     95  //  SetVerboseLevel(1);
    9296  corr = G4LossTableManager::Instance()->EmCorrections();
    9397}
     
    100104//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    101105
     106G4bool G4ionIonisation::IsApplicable(const G4ParticleDefinition& p)
     107{
     108  return (p.GetPDGCharge() != 0.0 && !p.IsShortLived() &&
     109          p.GetParticleType() == "nucleus");
     110}
     111
     112//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     113
     114G4double G4ionIonisation::MinPrimaryEnergy(const G4ParticleDefinition* p,
     115                                           const G4Material*,
     116                                           G4double cut)
     117{
     118  return
     119    p->GetPDGMass()*(std::sqrt(1. + 0.5*cut/CLHEP::electron_mass_c2) - 1.0);
     120}
     121
     122//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
     123
    102124void G4ionIonisation::InitialiseEnergyLossProcess(
    103125                      const G4ParticleDefinition* part,
    104126                      const G4ParticleDefinition* bpart)
    105127{
    106   if(isInitialised) return;
     128  const G4ParticleDefinition* ion = G4GenericIon::GenericIon();
    107129
    108   theParticle = part;
     130  if(!isInitialised) {
    109131
    110   if(part == bpart || part == G4GenericIon::GenericIon()) theBaseParticle = 0;
    111   else if(bpart == 0) theBaseParticle = G4GenericIon::GenericIon();
    112   else                theBaseParticle = bpart;
     132    theParticle = part;
     133    //G4String pname = part->GetParticleName();
    113134
    114   SetBaseParticle(theBaseParticle);
    115   SetSecondaryParticle(G4Electron::Electron());
     135    // define base particle
     136    const G4ParticleDefinition* theBaseParticle = 0;
    116137
    117   if(theBaseParticle) baseMass = theBaseParticle->GetPDGMass();
    118   else                baseMass = theParticle->GetPDGMass();
    119  
    120   if (!EmModel(1)) SetEmModel(new G4BraggIonModel(),1);
    121   EmModel(1)->SetLowEnergyLimit(100*eV);
    122   eth = 2.0*MeV; 
    123   EmModel(1)->SetHighEnergyLimit(eth);
    124   if (!FluctModel()) SetFluctModel(new G4IonFluctuations());
    125   AddEmModel(1, EmModel(1), FluctModel());
     138    if(part == ion)     theBaseParticle = 0;
     139    else if(bpart == 0) theBaseParticle = ion;
     140    else                theBaseParticle = bpart;
    126141
    127   if (!EmModel(2)) SetEmModel(new G4BetheBlochModel(),2); 
    128   EmModel(2)->SetLowEnergyLimit(eth);
    129   EmModel(2)->SetHighEnergyLimit(100*TeV);
    130   AddEmModel(2, EmModel(2), FluctModel());   
     142    SetBaseParticle(theBaseParticle);
     143    SetSecondaryParticle(G4Electron::Electron());
    131144
    132   // Add ion stoping tables for Generic Ion
    133   if(part == G4GenericIon::GenericIon()) {
    134     G4WaterStopping  ws(corr);
    135     effCharge = corr->GetIonEffectiveCharge(EmModel(1));
    136   } else {
    137     effCharge = corr->GetIonEffectiveCharge(0);
     145    if (!EmModel(1)) SetEmModel(new G4BraggIonModel(), 1);
     146    EmModel(1)->SetLowEnergyLimit(MinKinEnergy());
     147
     148    // model limit defined for protons
     149    eth = (EmModel(1)->HighEnergyLimit())*part->GetPDGMass()/proton_mass_c2;
     150    EmModel(1)->SetHighEnergyLimit(eth);
     151
     152    if (!FluctModel()) SetFluctModel(new G4IonFluctuations());
     153    AddEmModel(1, EmModel(1), FluctModel());
     154
     155    if (!EmModel(2)) SetEmModel(new G4BetheBlochModel(),2); 
     156    EmModel(2)->SetLowEnergyLimit(eth);
     157    EmModel(2)->SetHighEnergyLimit(MaxKinEnergy());
     158    AddEmModel(2, EmModel(2), FluctModel());   
     159
     160    // Add ion stoping tables for Generic Ion
     161    if(part == ion) {
     162      G4WaterStopping  ws(corr);
     163      corr->SetIonisationModels(EmModel(1),EmModel(2));
     164    }
     165    isInitialised = true;
    138166  }
    139 
    140   isInitialised = true;
     167  // reinitialisation of corrections for the new run
     168  EmModel(1)->ActivateNuclearStopping(nuclearStopping);
     169  EmModel(2)->ActivateNuclearStopping(nuclearStopping);
     170  if(part == ion) corr->InitialiseForNewRun();
    141171}
    142172
     
    145175void G4ionIonisation::PrintInfo()
    146176{
    147   if(EmModel(1) && EmModel(2))
    148     G4cout << "      Scaling relation is used from proton dE/dx and range."
    149            << "\n      Delta cross sections and sampling from "
    150            << EmModel(2)->GetName() << " model for scaled energy > "
    151            << eth/MeV << " MeV"
    152            << "\n      Parametrisation from "
    153            << EmModel(1)->GetName() << " for protons below."
    154            << " NuclearStopping= " << nuclearStopping
    155            << G4endl;   
    156   if (stopDataActive)
    157     G4cout << "\n      Stopping Power data for "
     177  if (stopDataActive && G4GenericIon::GenericIon() == theParticle) {
     178    G4cout << "      Stopping Power data for "
    158179           << corr->GetNumberOfStoppingVectors()
    159            << " ion/material pairs are used."
     180           << " ion/material pairs, nuclearStopping: " << nuclearStopping
    160181           << G4endl;
     182  }
    161183}
    162184
     
    165187void G4ionIonisation::AddStoppingData(G4int Z, G4int A,
    166188                                      const G4String& mname,
    167                                       G4PhysicsVector& dVector)
     189                                      G4PhysicsVector* dVector)
    168190{
    169191  corr->AddStoppingData(Z, A, mname, dVector);
     
    171193
    172194//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    173 
    174 void G4ionIonisation::CorrectionsAlongStep(const G4MaterialCutsCouple* couple,
    175                                            const G4DynamicParticle* dp,
    176                                            G4double& eloss,
    177                                            G4double& s)
    178 {
    179   const G4ParticleDefinition* part = dp->GetDefinition();
    180   const G4Material* mat = couple->GetMaterial();
    181   if(eloss < preKinEnergy) {
    182     //    G4cout << "e= " << preKinEnergy << " ratio= " << massRatio << " eth= " << eth<<  G4endl;
    183     if(preKinEnergy*massRatio > eth)
    184       eloss += s*corr->HighOrderCorrections(part,mat,preKinEnergy);
    185     else {
    186 
    187       if(stopDataActive)
    188         eloss *= corr->EffectiveChargeCorrection(part,mat,preKinEnergy);
    189 
    190     }
    191     fParticleChange.SetProposedCharge(effCharge->EffectiveCharge(part,
    192                                       mat,preKinEnergy-eloss));
    193   }
    194   if(nuclearStopping && preKinEnergy*massRatio < 50.*eth*charge2) {
    195     G4double nloss = s*corr->NuclearDEDX(part,mat,preKinEnergy - eloss*0.5);
    196     eloss += nloss;
    197     //  G4cout << "G4ionIonisation::CorrectionsAlongStep: e= " << preKinEnergy
    198     //     << " de= " << eloss << " NIEL= " << nloss << G4endl;
    199     fParticleChange.ProposeNonIonizingEnergyDeposit(nloss);
    200   }
    201 }
    202 
    203 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
Note: See TracChangeset for help on using the changeset viewer.