Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4DNABornExcitationModel2.cc
이 파일의 문서화 페이지로 가기
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: G4DNABornExcitationModel2.cc 90054 2015-05-11 18:47:32Z matkara $
27 //
28 
30 #include "G4SystemOfUnits.hh"
31 #include "G4DNAChemistryManager.hh"
33 #include "G4PhysicsTable.hh"
34 #include "G4PhysicsVector.hh"
35 #include "G4UnitsTable.hh"
36 #include <map>
37 
38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
39 
40 using namespace std;
41 
42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
43 
45  const G4String& nam) :
46  G4VEmModel(nam), isInitialised(false), fTableData(0)
47 {
49  fHighEnergy = 0;
50  fLowEnergy = 0;
52 
53  verboseLevel = 0;
54  // Verbosity scale:
55  // 0 = nothing
56  // 1 = warning for energy non-conservation
57  // 2 = details of energy budget
58  // 3 = calculation of cross sections, file openings, sampling of atoms
59  // 4 = entering in methods
60 
61  if (verboseLevel > 0)
62  {
63  G4cout << "Born excitation model is constructed " << G4endl;
64  }
67  fTotalXS = 0;
68  fTableData = 0;
69 
70  // Selection of stationary mode
71 
72  statCode = false;
73 }
74 
75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
76 
78 {
79  // Cross section
80  if (fTableData)
81  delete fTableData;
82 }
83 
84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
85 
87  const G4DataVector& /*cuts*/)
88 {
89 
90  if (verboseLevel > 3)
91  {
92  G4cout << "Calling G4DNABornExcitationModel2::Initialise()" << G4endl;
93  }
94 
95  if(fParticleDefinition != 0 && fParticleDefinition != particle)
96  {
97  G4Exception("G4DNABornExcitationModel2::Initialise","em0001",
98  FatalException,"Model already initialized for another particle type.");
99  }
100 
101  fParticleDefinition = particle;
102 
103  std::ostringstream fullFileName;
104  char *path = getenv("G4LEDATA");
105 
106  if(G4String(path) == "")
107  {
108  G4Exception("G4DNABornExcitationModel2::Initialise","G4LEDATA-CHECK",
109  FatalException, "G4LEDATA not defined in environment variables");
110  }
111 
112  fullFileName << path;
113 
114  if(particle->GetParticleName() == "e-")
115  {
116  fullFileName << "/dna/bornExcitation-e.dat";
117  fLowEnergy = 9*eV;
118  fHighEnergy = 1*MeV;
119  }
120  else if(particle->GetParticleName() == "proton")
121  {
122  fullFileName << "/dna/bornExcitation-p.dat";
123  fLowEnergy = 500. * keV;
124  fHighEnergy = 100. * MeV;
125  }
126 
129 
130 // G4double scaleFactor = (1.e-22 / 3.343) * m*m;
131 
132  fTableData = new G4PhysicsTable();
133  fTableData->RetrievePhysicsTable(fullFileName.str().c_str(), true);
134  for(size_t level = 0; level<fTableData->size(); ++level)
135  {
136 // (*fTableData)(level)->ScaleVector(1,scaleFactor);
137  (*fTableData)(level)->SetSpline(true);
138  }
139 
140  size_t finalBin_i = 2000;
141  G4double E_min = fLowEnergy;
142  G4double E_max = fHighEnergy;
143  fTotalXS = new G4PhysicsLogVector(E_min, E_max, finalBin_i);
144  fTotalXS->SetSpline(true);
146  G4double finalXS;
147 
148  for(size_t energy_i = 0; energy_i < finalBin_i; ++energy_i)
149  {
150  energy = fTotalXS->Energy(energy_i);
151  finalXS = 0;
152 
153  for(size_t level = 0; level<fTableData->size(); ++level)
154  {
155  finalXS += (*fTableData)(level)->Value(energy);
156  }
157  fTotalXS->PutValue(energy_i, finalXS);
158 // G4cout << "energy = " << energy << " " << fTotalXS->Value(energy)
159 // << " " << energy_i << " " << finalXS << G4endl;
160  }
161 
162 // for(energy = LowEnergyLimit() ; energy < HighEnergyLimit() ; energy += 1*pow(10,log10(energy)))
163 // {
164 // G4cout << "energy = " << energy << " " << fTotalXS->Value(energy) << G4endl;
165 // }
166 
167  if( verboseLevel>0 )
168  {
169  G4cout << "Born excitation model is initialized " << G4endl
170  << "Energy range: "
171  << LowEnergyLimit() / eV << " eV - "
172  << HighEnergyLimit() / keV << " keV for "
173  << particle->GetParticleName()
174  << G4endl;
175  }
176 
177  // Initialize water density pointer
179 
180  if (isInitialised)
181  { return;}
183  isInitialised = true;
184 }
185 
186 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
187 
189  const G4ParticleDefinition* particleDefinition,
190  G4double ekin,
191  G4double,
192  G4double)
193 {
194  if (verboseLevel > 3)
195  {
196  G4cout << "Calling CrossSectionPerVolume() of G4DNABornExcitationModel2"
197  << G4endl;
198  }
199 
200  if(particleDefinition != fParticleDefinition) return 0;
201 
202  // Calculate total cross section for model
203 
204  G4double sigma=0;
205 
206  G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
207 
208  if(waterDensity!= 0.0)
209  {
210  if (ekin >= fLowEnergy && ekin < fHighEnergy)
211  {
212  sigma = fTotalXS->Value(ekin, fLastBinCallForFinalXS);
213 
214 // for(size_t i = 0; i < 5; ++i)
215 // sigma += (*fTableData)[i]->Value(ekin);
216 
217  if(sigma == 0)
218  {
219  G4cerr << "PROBLEM SIGMA = 0 at " << G4BestUnit(ekin, "Energy")<< G4endl;
220  }
221  }
222 
223  if (verboseLevel > 2)
224  {
225  G4cout << "__________________________________" << G4endl;
226  G4cout << "G4DNABornExcitationModel2 - XS INFO START" << G4endl;
227  G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
228  G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
229  G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
230  G4cout << "G4DNABornExcitationModel2 - XS INFO END" << G4endl;
231  }
232  } // if (waterMaterial)
233  else
234  {
235  G4cerr << "PROBLEM NO WATER MATERIAL ?" << G4endl;
236  }
237 
238  return sigma*waterDensity;
239 }
240 
241 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
242 
243 void G4DNABornExcitationModel2::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
244  const G4MaterialCutsCouple* /*couple*/,
245  const G4DynamicParticle* aDynamicParticle,
246  G4double,
247  G4double)
248 {
249 
250  if (verboseLevel > 3)
251  {
252  G4cout << "Calling SampleSecondaries() of G4DNABornExcitationModel2"
253  << G4endl;
254  }
255 
256  G4double k = aDynamicParticle->GetKineticEnergy();
257 
258  G4int level = RandomSelect(k);
259  G4double excitationEnergy = waterStructure.ExcitationEnergy(level);
260  G4double newEnergy = k - excitationEnergy;
261 
262  if (newEnergy > 0)
263  {
265 
268 
270  }
271 
272  const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
274  level,
275  theIncomingTrack);
276 }
277 
278 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
279 
281  G4int level,
282  const G4ParticleDefinition* particle,
284 {
285  if (fParticleDefinition != particle)
286  {
287  G4Exception("G4DNABornExcitationModel2::GetPartialCrossSection",
288  "bornParticleType",
290  "Model initialized for another particle type.");
291  }
292 
293  return (*fTableData)(level)->Value(kineticEnergy);
294 }
295 
296 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
297 
299 {
300  const size_t n(fTableData->size());
301  size_t i(n);
302 
304 
305  value *= G4UniformRand();
306  i = n;
307 
308  G4double partialXS;
309 
310  while (i > 0)
311  {
312  i--;
313 
314  partialXS = (*fTableData)(i)->Value(k);
315  if (partialXS > value)
316  {
317  return i;
318  }
319  value -= partialXS;
320  }
321 
322  return 0;
323 }
324 
G4double Energy(size_t index) const
static constexpr double MeV
Definition: G4SIunits.hh:214
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:707
static constexpr double keV
Definition: G4SIunits.hh:216
size_t GetIndex() const
Definition: G4Material.hh:262
#define G4endl
Definition: G4ios.hh:61
const G4ThreeVector & GetMomentumDirection() const
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
const G4String & GetParticleName() const
G4DNABornExcitationModel2(const G4ParticleDefinition *p=0, const G4String &nam="DNABornExcitationModel")
void SetProposedKineticEnergy(G4double proposedKinEnergy)
G4double Value(G4double theEnergy, size_t &lastidx) const
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:608
const std::vector< G4double > * fpMolWaterDensity
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:616
const std::vector< double > * GetNumMolPerVolTableFor(const G4Material *) const
Retrieve a table of molecular densities (number of molecules per unit volume) in the G4 unit system f...
G4DNAWaterExcitationStructure waterStructure
double G4double
Definition: G4Types.hh:76
const XML_Char int const XML_Char * value
Definition: expat.h:331
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &=*(new G4DataVector()))
static G4DNAMolecularMaterial * Instance()
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
#define G4UniformRand()
Definition: Randomize.hh:53
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
void SetSpline(G4bool)
double energy
Definition: plottest35.C:25
static constexpr double eV
Definition: G4SIunits.hh:215
const G4Track * GetCurrentTrack() const
G4GLOB_DLL std::ostream G4cerr
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
G4bool RetrievePhysicsTable(const G4String &filename, G4bool ascii=false)
int G4int
Definition: G4Types.hh:78
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:714
G4ParticleChangeForGamma * fParticleChangeForGamma
G4double GetKineticEnergy() const
static constexpr double cm
Definition: G4SIunits.hh:119
G4GLOB_DLL std::ostream G4cout
Char_t n[5]
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:131
const G4ParticleDefinition * fParticleDefinition
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:609
virtual G4double GetPartialCrossSection(const G4Material *, G4int level, const G4ParticleDefinition *, G4double kineticEnergy)
static G4DNAChemistryManager * Instance()
virtual G4double Value(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
Definition: G4VEmModel.cc:413
void PutValue(size_t index, G4double theValue)