source: trunk/source/processes/electromagnetic/lowenergy/src/G4DNABornExcitationModel.cc @ 1317

Last change on this file since 1317 was 1315, checked in by garnier, 14 years ago

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 9.7 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer                                           *
4// *                                                                  *
5// * The  Geant4 software  is  copyright of the Copyright Holders  of *
6// * the Geant4 Collaboration.  It is provided  under  the terms  and *
7// * conditions of the Geant4 Software License,  included in the file *
8// * LICENSE and available at  http://cern.ch/geant4/license .  These *
9// * include a list of copyright holders.                             *
10// *                                                                  *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work  make  any representation or  warranty, express or implied, *
14// * regarding  this  software system or assume any liability for its *
15// * use.  Please see the license in the file  LICENSE  and URL above *
16// * for the full disclaimer and the limitation of liability.         *
17// *                                                                  *
18// * This  code  implementation is the result of  the  scientific and *
19// * technical work of the GEANT4 collaboration.                      *
20// * By using,  copying,  modifying or  distributing the software (or *
21// * any work based  on the software)  you  agree  to acknowledge its *
22// * use  in  resulting  scientific  publications,  and indicate your *
23// * acceptance of all terms of the Geant4 Software license.          *
24// ********************************************************************
25//
26// $Id: G4DNABornExcitationModel.cc,v 1.9 2010/03/27 11:32:41 sincerti Exp $
27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
28//
29
30#include "G4DNABornExcitationModel.hh"
31
32//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
33
34using namespace std;
35
36//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
37
38G4DNABornExcitationModel::G4DNABornExcitationModel(const G4ParticleDefinition*,
39                                             const G4String& nam)
40:G4VEmModel(nam),isInitialised(false)
41{
42  verboseLevel= 0;
43  // Verbosity scale:
44  // 0 = nothing
45  // 1 = warning for energy non-conservation
46  // 2 = details of energy budget
47  // 3 = calculation of cross sections, file openings, sampling of atoms
48  // 4 = entering in methods
49 
50  if( verboseLevel>0 ) 
51  { 
52    G4cout << "Born excitation model is constructed " << G4endl;
53  }
54 
55}
56
57//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
58
59G4DNABornExcitationModel::~G4DNABornExcitationModel()
60{ 
61  // Cross section
62 
63  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
64  for (pos = tableData.begin(); pos != tableData.end(); ++pos)
65  {
66    G4DNACrossSectionDataSet* table = pos->second;
67    delete table;
68  }
69
70}
71
72//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
73
74void G4DNABornExcitationModel::Initialise(const G4ParticleDefinition* particle,
75                                       const G4DataVector& /*cuts*/)
76{
77
78  if (verboseLevel > 3)
79    G4cout << "Calling G4DNABornExcitationModel::Initialise()" << G4endl;
80
81  // Energy limits
82  // Energy limits
83
84  G4String fileElectron("dna/sigma_excitation_e_born");
85  G4String fileProton("dna/sigma_excitation_p_born");
86
87  G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
88  G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
89
90  G4String electron;
91  G4String proton;
92 
93  G4double scaleFactor = (1.e-22 / 3.343) * m*m;
94
95  if (electronDef != 0)
96  {
97    electron = electronDef->GetParticleName();
98
99    tableFile[electron] = fileElectron;
100
101    lowEnergyLimit[electron] = 9. * eV; 
102    highEnergyLimit[electron] = 1. * MeV;
103
104    // Cross section
105   
106    G4DNACrossSectionDataSet* tableE = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
107    tableE->LoadData(fileElectron);
108     
109    tableData[electron] = tableE;
110   
111  }
112  else
113  {
114    G4Exception("G4DNABornExcitationModel::Initialise(): electron is not defined");
115  }
116
117  if (protonDef != 0)
118  {
119    proton = protonDef->GetParticleName();
120
121    tableFile[proton] = fileProton;
122
123    lowEnergyLimit[proton] = 500. * keV;
124    highEnergyLimit[proton] = 100. * MeV;
125
126    // Cross section
127   
128    G4DNACrossSectionDataSet* tableP = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
129    tableP->LoadData(fileProton);
130     
131    tableData[proton] = tableP;
132     
133  }
134  else
135  {
136    G4Exception("G4DNABornExcitationModel::Initialise(): proton is not defined");
137  }
138
139  if (particle==electronDef) 
140  {
141    SetLowEnergyLimit(lowEnergyLimit[electron]);
142    SetHighEnergyLimit(highEnergyLimit[electron]);
143  }
144
145  if (particle==protonDef) 
146  {
147    SetLowEnergyLimit(lowEnergyLimit[proton]);
148    SetHighEnergyLimit(highEnergyLimit[proton]);
149  }
150
151  if( verboseLevel>0 ) 
152  { 
153    G4cout << "Born excitation model is initialized " << G4endl
154           << "Energy range: "
155           << LowEnergyLimit() / eV << " eV - "
156           << HighEnergyLimit() / keV << " keV for "
157           << particle->GetParticleName()
158           << G4endl;
159  }
160 
161 
162  if(!isInitialised) 
163  {
164    isInitialised = true;
165 
166    if(pParticleChange)
167      fParticleChangeForGamma = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
168    else
169      fParticleChangeForGamma = new G4ParticleChangeForGamma();
170  }   
171
172  // InitialiseElementSelectors(particle,cuts);
173
174}
175
176//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
177
178G4double G4DNABornExcitationModel::CrossSectionPerVolume(const G4Material* material,
179                                           const G4ParticleDefinition* particleDefinition,
180                                           G4double ekin,
181                                           G4double,
182                                           G4double)
183{
184  if (verboseLevel > 3)
185    G4cout << "Calling CrossSectionPerVolume() of G4DNABornExcitationModel" << G4endl;
186
187  if (
188      particleDefinition != G4Proton::ProtonDefinition()
189      &&
190      particleDefinition != G4Electron::ElectronDefinition()
191     )
192           
193    return 0;
194
195 // Calculate total cross section for model
196
197  G4double lowLim = 0;
198  G4double highLim = 0;
199  G4double sigma=0;
200
201  if (material->GetName() == "G4_WATER")
202  {
203    const G4String& particleName = particleDefinition->GetParticleName();
204 
205    std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
206    pos1 = lowEnergyLimit.find(particleName);
207    if (pos1 != lowEnergyLimit.end())
208    {
209      lowLim = pos1->second;
210    }
211 
212    std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
213    pos2 = highEnergyLimit.find(particleName);
214    if (pos2 != highEnergyLimit.end())
215    {
216      highLim = pos2->second;
217    }
218
219    if (ekin >= lowLim && ekin < highLim)
220    {
221      std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
222      pos = tableData.find(particleName);
223       
224      if (pos != tableData.end())
225      {
226        G4DNACrossSectionDataSet* table = pos->second;
227        if (table != 0)
228        {
229          sigma = table->FindValue(ekin);
230        }
231      }
232      else
233      {
234        G4Exception("G4DNABornExcitationModel::CrossSectionPerVolume: attempting to calculate cross section for wrong particle");
235      }
236    }
237
238    if (verboseLevel > 3)
239    {
240      G4cout << "---> Kinetic energy(eV)=" << ekin/eV << G4endl;
241      G4cout << " - Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
242      G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
243    } 
244 
245  } // if (waterMaterial)
246 
247 return sigma*material->GetAtomicNumDensityVector()[1];           
248           
249
250}
251
252//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
253
254void G4DNABornExcitationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
255                                              const G4MaterialCutsCouple* /*couple*/,
256                                              const G4DynamicParticle* aDynamicParticle,
257                                              G4double,
258                                              G4double)
259{
260
261  if (verboseLevel > 3)
262    G4cout << "Calling SampleSecondaries() of G4DNABornExcitationModel" << G4endl;
263
264  G4double k = aDynamicParticle->GetKineticEnergy();
265 
266  const G4String& particleName = aDynamicParticle->GetDefinition()->GetParticleName();
267
268  G4int level = RandomSelect(k,particleName);
269  G4double excitationEnergy = waterStructure.ExcitationEnergy(level);
270  G4double newEnergy = k - excitationEnergy;
271 
272  if (newEnergy > 0)
273  {
274      fParticleChangeForGamma->ProposeMomentumDirection(aDynamicParticle->GetMomentumDirection());
275      fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
276      fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
277  }
278
279}
280
281//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
282
283G4int G4DNABornExcitationModel::RandomSelect(G4double k, const G4String& particle)
284{   
285  G4int level = 0;
286
287  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
288  pos = tableData.find(particle);
289
290  if (pos != tableData.end())
291  {
292    G4DNACrossSectionDataSet* table = pos->second;
293
294    if (table != 0)
295    {
296      G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
297      const size_t n(table->NumberOfComponents());
298      size_t i(n);
299      G4double value = 0.;
300           
301      while (i>0)
302      { 
303        i--;
304        valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
305        value += valuesBuffer[i];
306      }
307           
308      value *= G4UniformRand();
309   
310      i = n;
311           
312      while (i > 0)
313      {
314        i--;
315
316        if (valuesBuffer[i] > value)
317        {
318          delete[] valuesBuffer;
319          return i;
320        }
321        value -= valuesBuffer[i];
322      }
323           
324      if (valuesBuffer) delete[] valuesBuffer;
325   
326    }
327  }
328  else
329  {
330    G4Exception("G4DNABornExcitationModel::RandomSelect attempting to calculate cross section for wrong particle");
331  }
332  return level;
333}
334
335
Note: See TracBrowser for help on using the repository browser.