Message: Ionization Energy Not Logged In (login)
 Next-in-Thread Next-in-Thread
 Next-in-Forum Next-in-Forum

None Ionization Energy 

Forum: Event and Track Management
Date: 12 Feb, 2012
From: tbjc1magic <tbjc1magic>

Hi,

 I am a new programmer.Currently, I am trying to simulate how much ionization energy particles will deposit when they go through different gas.

 1.Therefore, I read Geant4 example "TestEm18", which accumulates all charged secondary particles' kinetic energy. It uses this energy as ionization energy.

 2.I also tried another method: for each step, g4step has the information GetTotalEnergyDeposit() and GetNonIonizingEnergyDeposit(); I take their difference, and add up those difference from each step. I think this can be ionization energy.

 However, the results told me the first method is smaller than the second one, why is that and which method is correct?

 Attachment is the result. eLoss/primary is the ionization energy from 1st method; ionizationenergy is the one from 2ed method.

 By the way, can anybody tell me what "TestEm18" tries to calculate in the following code:

void RunAction::BeginOfRunAction(const G4Run* run)
{
  G4cout << "### Run " << run->GetRunID() << " start." << G4endl;

  //initialisation
  //
  energyDeposit = 0.;

  nbCharged = nbNeutral = 0;
  energyCharged = energyNeutral = 0.;
	ionizationenergy=0;
  emin[0] = emin[1] = DBL_MAX;
  emax[0] = emax[1] = 0.;    

  nbSteps = 0;
  trackLength = 0.; 

  histoManager->book();

  // do not save Rndm status
  G4RunManager::GetRunManager()->SetRandomNumberStore(false);
  CLHEP::HepRandom::showEngineStatus();
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void RunAction::EndOfRunAction(const G4Run* aRun)
{
  G4int nbEvents = aRun->GetNumberOfEvent();
  if (nbEvents == 0) return;

  G4Material* material = detector->GetMaterial();
  G4double length  = detector->GetSize();
  G4double density = material->GetDensity();

  G4ParticleDefinition* particle = primary->GetParticleGun()
                                          ->GetParticleDefinition();
  G4String partName = particle->GetParticleName();
  G4double eprimary = primary->GetParticleGun()->GetParticleEnergy();

  G4int prec = G4cout.precision(3);
  G4cout << "\n ======================== run summary ======================\n";  
  G4cout << "\n The run was " << nbEvents << " " << partName << " of "
         << G4BestUnit(eprimary,"Energy") << " through " 
	 << G4BestUnit(length,"Length") << " of "
	 << material->GetName() << " (density: " 
	 << G4BestUnit(density,"Volumic Mass") << ")";
  G4cout << "\n ===========================================================\n";
  G4cout << G4endl;

  histoManager->save();

  if (particle->GetPDGCharge() == 0.) return;

  G4cout.precision(5);

  //track length
  //
  G4double trackLPerEvent = trackLength/nbEvents;
  G4double nbStepPerEvent = double(nbSteps)/nbEvents;
  G4double stepSize = trackLength/nbSteps;

  G4cout 
    << "\n trackLength= " 
    << G4BestUnit(trackLPerEvent, "Length")
    << "\t nb of steps= " << nbStepPerEvent
    << "  stepSize= " << G4BestUnit(stepSize, "Length")
    << G4endl;

  //charged secondaries (ionization, direct pair production)
  //
  G4double energyPerEvent = energyCharged/nbEvents;
  G4double nbPerEvent = double(nbCharged)/nbEvents;
  G4double meanEkin = 0.;
  if (nbCharged) meanEkin = energyCharged/nbCharged;

  G4cout 
    << "\n d-rays  : eLoss/primary= " 
    << G4BestUnit(energyPerEvent, "Energy")
    << "\t  nb of d-rays= " << nbPerEvent
    << "  <Tkin>= " << G4BestUnit(meanEkin, "Energy")
    << "  Tmin= "   << G4BestUnit(emin[0],  "Energy")
    << "  Tmax= "   << G4BestUnit(emax[0],  "Energy")
    << G4endl;

// G4cout<<"ionizationenergy is "<<G4BestUnit(ionizationenergy/nbEvents, "Energy")<<G4endl;

  //neutral secondaries (bremsstrahlung, pixe)
  //
  energyPerEvent = energyNeutral/nbEvents;
  nbPerEvent = double(nbNeutral)/nbEvents;
  meanEkin = 0.;
  if (nbNeutral) meanEkin = energyNeutral/nbNeutral;

  G4cout 
    << "\n gamma   : eLoss/primary= " 
    << G4BestUnit(energyPerEvent, "Energy")
    << "\t  nb of gammas= " << nbPerEvent
    << "  <Tkin>= " << G4BestUnit(meanEkin, "Energy")
    << "  Tmin= "   << G4BestUnit(emin[1],  "Energy")
    << "  Tmax= "   << G4BestUnit(emax[1],  "Energy")
    << G4endl;

  G4EmCalculator emCal;

  //local energy deposit
  //
  energyPerEvent = energyDeposit/nbEvents;
  //
  G4double r0  = emCal.GetRangeFromRestricteDEDX(eprimary,particle,material);  
  G4double r1 = r0 - trackLPerEvent;
  G4double etry = eprimary - energyPerEvent;  
  G4double efinal = 0.;
  if (r1 > 0.) efinal = GetEnergyFromRestrictedRange(r1,particle,material,etry);
  G4double dEtable = eprimary - efinal;
  G4double ratio = 0.;
  if (dEtable > 0.) ratio = energyPerEvent/dEtable;

  G4cout 
    << "\n deposit : eLoss/primary= " 
    << G4BestUnit(energyPerEvent, "Energy")
    << "\t <dEcut > table= " 
    << G4BestUnit(dEtable, "Energy")
    << "   ---> simul/reference= " << ratio           
    << G4endl;

  //total energy transferred
  //
  G4double energyTotal = energyDeposit + energyCharged + energyNeutral;
  energyPerEvent = energyTotal/nbEvents;
  //
  r0  = emCal.GetCSDARange(eprimary,particle,material);  
  r1 = r0 - trackLPerEvent;
  etry = eprimary - energyPerEvent;
  efinal = 0.;
  if (r1 > 0.) efinal = GetEnergyFromCSDARange(r1,particle,material,etry);
  dEtable = eprimary - efinal;
  ratio = 0.;
  if (dEtable > 0.) ratio = energyPerEvent/dEtable;

  G4cout 
    << "\n total   : eLoss/primary= " 
    << G4BestUnit(energyPerEvent, "Energy")
    << "\t <dEfull> table= " 
    << G4BestUnit(dEtable, "Energy")
    << "   ---> simul/reference= " << ratio           
    << G4endl; 

  G4cout.precision(prec);

  // show Rndm status
  CLHEP::HepRandom::showEngineStatus();
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

G4double RunAction::GetEnergyFromRestrictedRange(G4double range,
            G4ParticleDefinition* particle, G4Material* material, G4double Etry)
{
  G4EmCalculator emCal;

  G4double Energy = Etry, dE = 0., dEdx;
  G4double r, dr;
  G4double err  = 1., errmax = 0.00001;
  G4int    iter = 0 , itermax = 10;  
  while (err > errmax && iter < itermax) {
    iter++;
    Energy -= dE;
    r = emCal.GetRangeFromRestricteDEDX(Energy,particle,material);
    dr = r - range;          
    dEdx = emCal.GetDEDX(Energy,particle,material);
    dE = dEdx*dr;
    err = std::abs(dE)/Energy;    
  }
  if (iter == itermax) {
    G4cout 
    << "\n  ---> warning: RunAction::GetEnergyFromRestRange() did not converge"
    << "   Etry = " << G4BestUnit(Etry,"Energy")
    << "   Energy = " << G4BestUnit(Energy,"Energy")
    << "   err = " << err
    << "   iter = " << iter << G4endl;
  }	 

  return Energy;	 
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

G4double RunAction::GetEnergyFromCSDARange(G4double range,
            G4ParticleDefinition* particle, G4Material* material, G4double Etry)
{
  G4EmCalculator emCal;

  G4double Energy = Etry, dE = 0., dEdx;
  G4double r, dr;
  G4double err  = 1., errmax = 0.00001;
  G4int    iter = 0 , itermax = 10;  
  while (err > errmax && iter < itermax) {
    iter++;
    Energy -= dE;
    r = emCal.GetCSDARange(Energy,particle,material);
    dr = r - range;          
    dEdx = emCal.ComputeTotalDEDX(Energy,particle,material);
    dE = dEdx*dr;
    err = std::abs(dE)/Energy;
  }
  if (iter == itermax) {
    G4cout 
    << "\n  ---> warning: RunAction::GetEnergyFromCSDARange() did not converge"
    << "   Etry = " << G4BestUnit(Etry,"Energy")
    << "   Energy = " << G4BestUnit(Energy,"Energy")
    << "   err = " << err
    << "   iter = " << iter << G4endl;
  }	 

  return Energy;	 
}

   Attachment:
      http://hypernews.slac.stanford.edu/HyperNews/geant4/get/AUX/2012/02/12/14.18-91084-9B2-9817-5F9212FE763B.png

Inline Depth:
 1 1
 All All
Outline Depth:
 1 1
 2 2
 All All
Add message: (add)

1 Question: Re: Ionization Energy   (tbjc1magic - 13 Feb, 2012)
(_ Feedback: Re: Ionization Energy   (Gumplinger Peter - 13 Feb, 2012)
(_ Question: Re: Ionization Energy   (tbjc1magic - 13 Feb, 2012)
(_ Feedback: Re: Ionization Energy   (Gumplinger Peter - 14 Feb, 2012)
(_ Question: Re: Ionization Energy   (tbjc1magic - 14 Feb, 2012)
(_ Feedback: Re: Ionization Energy   (Gumplinger Peter - 14 Feb, 2012)
 Add Message Add Message
to: "Ionization Energy"

 Subscribe Subscribe

This site runs SLAC HyperNews version 1.11-slac-98, derived from the original HyperNews


[ Geant 4 Home | Geant 4 HyperNews | Search | Request New Forum | Feedback ]