Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/dna/models/src/G4DNAQuinnPlasmonExcitationModel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/electromagnetic/dna/models/src/G4DNAQuinnPlasmonExcitationModel.cc (Version 11.3.0) and /processes/electromagnetic/dna/models/src/G4DNAQuinnPlasmonExcitationModel.cc (Version 11.0)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 // Created on 2016/04/08                           26 // Created on 2016/04/08
 27 //                                                 27 //
 28 // Authors: D. Sakata, S. Incerti                  28 // Authors: D. Sakata, S. Incerti
 29 //                                                 29 //
 30 // This class perform transmission term of vol     30 // This class perform transmission term of volume plasmon excitation,
 31 // based on Quinn Model, see Phys. Rev. vol 12     31 // based on Quinn Model, see Phys. Rev. vol 126, number 4 (1962)
 32                                                    32 
 33 #include "G4DNAQuinnPlasmonExcitationModel.hh"     33 #include "G4DNAQuinnPlasmonExcitationModel.hh"
 34 #include "G4SystemOfUnits.hh"                      34 #include "G4SystemOfUnits.hh"
 35 #include "G4RandomDirection.hh"                    35 #include "G4RandomDirection.hh"
 36                                                    36 
 37 //....oooOO0OOooo........oooOO0OOooo........oo     37 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 38                                                    38 
 39 using namespace std;                               39 using namespace std;
 40                                                    40 
 41 //....oooOO0OOooo........oooOO0OOooo........oo     41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 42                                                    42 
 43 G4DNAQuinnPlasmonExcitationModel::G4DNAQuinnPl     43 G4DNAQuinnPlasmonExcitationModel::G4DNAQuinnPlasmonExcitationModel
 44                         (const G4ParticleDefin     44                         (const G4ParticleDefinition*,
 45                          const G4String& nam):     45                          const G4String& nam):
 46 G4VEmModel(nam)                                <<  46 G4VEmModel(nam), isInitialised(false)
 47 {                                                  47 {
 48   fpMaterialDensity       = nullptr;           <<  48   fpMaterialDensity       = 0;
 49   fLowEnergyLimit         =  10  *  eV;            49   fLowEnergyLimit         =  10  *  eV;
 50   fHighEnergyLimit        =  1.0 * GeV;            50   fHighEnergyLimit        =  1.0 * GeV;
 51                                                    51 
 52   for(int & i : nValenceElectron) i=0;         <<  52   for(G4int i=0;i<100;i++) nValenceElectron[i]=0;
 53                                                    53 
 54   verboseLevel = 0;                                54   verboseLevel = 0;
 55                                                    55 
 56   if (verboseLevel > 0)                            56   if (verboseLevel > 0)
 57   {                                                57   {
 58     G4cout << "Quinn plasmon excitation model      58     G4cout << "Quinn plasmon excitation model is constructed " << G4endl;
 59   }                                                59   }
 60   fParticleChangeForGamma = nullptr;           <<  60   fParticleChangeForGamma = 0;
 61   statCode                = false;                 61   statCode                = false;
 62 }                                                  62 }
 63                                                    63 
 64 //....oooOO0OOooo........oooOO0OOooo........oo     64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 65                                                    65 
 66 G4DNAQuinnPlasmonExcitationModel::~G4DNAQuinnP     66 G4DNAQuinnPlasmonExcitationModel::~G4DNAQuinnPlasmonExcitationModel()
 67 = default;                                     <<  67 {
                                                   >>  68 }
 68                                                    69 
 69 //....oooOO0OOooo........oooOO0OOooo........oo     70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 70                                                    71 
 71 void G4DNAQuinnPlasmonExcitationModel::Initial     72 void G4DNAQuinnPlasmonExcitationModel::Initialise
 72                         (const G4ParticleDefin     73                         (const G4ParticleDefinition* particle,
 73                          const G4DataVector& /     74                          const G4DataVector& /*cuts*/)
 74 {                                                  75 {
 75   for(int & i : nValenceElectron) i=0;         <<  76   for(G4int i=0;i<100;i++) nValenceElectron[i]=0;
 76                                                    77 
 77   if (verboseLevel > 3)                            78   if (verboseLevel > 3)
 78   {                                                79   {
 79     G4cout <<                                      80     G4cout << 
 80            "Calling G4DNAQuinnPlasmonExcitatio     81            "Calling G4DNAQuinnPlasmonExcitationModel::Initialise()" 
 81            << G4endl;                              82            << G4endl;
 82   }                                                83   }
 83                                                    84 
 84   if(particle == G4Electron::ElectronDefinitio     85   if(particle == G4Electron::ElectronDefinition())
 85   {                                                86   {
 86     fLowEnergyLimit           =  10  *  eV;        87     fLowEnergyLimit           =  10  *  eV;
 87     fHighEnergyLimit          =  1.0 * GeV;        88     fHighEnergyLimit          =  1.0 * GeV;
 88   }                                                89   }
 89   else                                             90   else
 90   {                                                91   { 
 91     G4Exception("G4DNAQuinnPlasmonExcitationMo     92     G4Exception("G4DNAQuinnPlasmonExcitationModel::Initialise","em0001",
 92     FatalException,"Not defined for other part     93     FatalException,"Not defined for other particles than electrons.");
 93    return;                                         94    return;
 94   }                                                95   }
 95                                                    96 
 96   // Get Number of valence electrons               97   // Get Number of valence electrons
 97   G4ProductionCutsTable* theCoupleTable =          98   G4ProductionCutsTable* theCoupleTable = 
 98       G4ProductionCutsTable::GetProductionCuts     99       G4ProductionCutsTable::GetProductionCutsTable();
 99                                                   100       
100   auto  numOfCouples = (G4int)theCoupleTable-> << 101   G4int numOfCouples = theCoupleTable->GetTableSize();
101                                                   102   
102   for(G4int i=0;i<numOfCouples;i++){              103   for(G4int i=0;i<numOfCouples;i++){
103                                                   104     
104     const G4MaterialCutsCouple* couple =          105     const G4MaterialCutsCouple* couple = 
105           theCoupleTable->GetMaterialCutsCoupl    106           theCoupleTable->GetMaterialCutsCouple(i);
106                                                   107     
107     const G4Material* material = couple->GetMa    108     const G4Material* material = couple->GetMaterial();
108                                                   109     
109     const G4ElementVector* theElementVector =m    110     const G4ElementVector* theElementVector =material->GetElementVector();
110                                                   111     
111     std::size_t nelm = material->GetNumberOfEl << 112     G4int nelm = material->GetNumberOfElements();
112     if (nelm==1) // Protection: only for singl << 113     if (nelm==1){// Protection: only for single element
113     {                                          << 
114       G4int z = G4lrint((*theElementVector)[0]    114       G4int z = G4lrint((*theElementVector)[0]->GetZ());
115       if(z<=100)                               << 115       if(z<=100){nValenceElectron[z]  = GetNValenceElectron(z);}
116       {                                        << 
117         nValenceElectron[z]  = GetNValenceElec << 
118       }                                        << 
119       else                                     << 
120       {                                        << 
121         G4Exception("G4DNAQuinnPlasmonExcitati << 
122           FatalException,"The model is not app << 
123       }                                        << 
124     }                                             116     }
125     //for(G4int j=0;j<nelm;j++){                  117     //for(G4int j=0;j<nelm;j++){
126     //    G4int z=G4lrint((*theElementVector)[    118     //    G4int z=G4lrint((*theElementVector)[j]->GetZ());
127     //    if(z<=100){nValenceElectron[z]  = Ge    119     //    if(z<=100){nValenceElectron[z]  = GetNValenceElectron(z);}
128     //}                                           120     //}
129   }                                               121   }
130                                                   122 
131   if( verboseLevel>0 )                            123   if( verboseLevel>0 )
132   {                                               124   {
133     G4cout << "Quinn plasmon excitation model     125     G4cout << "Quinn plasmon excitation model is initialized " << G4endl
134     << "Energy range: "                           126     << "Energy range: "
135     << LowEnergyLimit() / eV << " eV - "          127     << LowEnergyLimit() / eV << " eV - "
136     << HighEnergyLimit() / keV << " keV for "     128     << HighEnergyLimit() / keV << " keV for "
137     << particle->GetParticleName()                129     << particle->GetParticleName()
138     << G4endl;                                    130     << G4endl;
139   }                                               131   }
140                                                   132 
141   if (isInitialised){return;}                     133   if (isInitialised){return;}
142   fParticleChangeForGamma = GetParticleChangeF    134   fParticleChangeForGamma = GetParticleChangeForGamma();
143   isInitialised = true;                           135   isInitialised = true;
144 }                                                 136 }
145                                                   137 
146 //....oooOO0OOooo........oooOO0OOooo........oo    138 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
147                                                   139 
148 G4double G4DNAQuinnPlasmonExcitationModel::Cro    140 G4double G4DNAQuinnPlasmonExcitationModel::CrossSectionPerVolume
149                          (const G4Material* ma    141                          (const G4Material* material,
150                           const G4ParticleDefi    142                           const G4ParticleDefinition* particleDefinition,
151                           G4double ekin,          143                           G4double ekin,
152                           G4double,               144                           G4double,
153                           G4double)               145                           G4double)
154 {                                                 146 {
155   if (verboseLevel > 3)                           147   if (verboseLevel > 3)
156   {                                               148   {
157     G4cout <<                                     149     G4cout << 
158          "Calling CrossSectionPerVolume() of G    150          "Calling CrossSectionPerVolume() of G4DNAQuinnPlasmonExcitationModel"
159            << G4endl;                             151            << G4endl;
160   }                                               152   }
161                                                   153 
162   // Protection: only for single element          154   // Protection: only for single element
163   if(material->GetNumberOfElements()>1) return    155   if(material->GetNumberOfElements()>1) return 0.;
164   G4double z = material->GetZ();                  156   G4double z = material->GetZ();
165                                                   157 
166   // Protection: only for Gold                    158   // Protection: only for Gold
167   if (z!=79){return 0.;}                          159   if (z!=79){return 0.;}
168                                                   160   
169                                                   161 
170   G4double sigma           = 0;                   162   G4double sigma           = 0;
171   G4double atomicNDensity = material->GetAtomi    163   G4double atomicNDensity = material->GetAtomicNumDensityVector()[0];
172                                                   164 
173   if(atomicNDensity!= 0.0)                        165   if(atomicNDensity!= 0.0)
174   {                                               166   {
175     if (ekin >= fLowEnergyLimit && ekin < fHig    167     if (ekin >= fLowEnergyLimit && ekin < fHighEnergyLimit)
176     {                                             168     {
177       sigma = GetCrossSection(material,particl    169       sigma = GetCrossSection(material,particleDefinition,ekin);
178     }                                             170     }
179                                                   171 
180     if (verboseLevel > 2)                         172     if (verboseLevel > 2)
181     {                                             173     {
182       G4cout<<"_______________________________    174       G4cout<<"__________________________________" << G4endl;
183       G4cout<<"=== G4DNAQuinnPlasmonExcitation    175       G4cout<<"=== G4DNAQuinnPlasmonExcitationModel - XS INFO START"<<G4endl;
184       G4cout<<"=== Kinetic energy (eV)=" << ek    176       G4cout<<"=== Kinetic energy (eV)=" << ekin/eV << " particle : " 
185             <<particleDefinition->GetParticleN    177             <<particleDefinition->GetParticleName() << G4endl;
186       G4cout<<"=== Cross section per atom for     178       G4cout<<"=== Cross section per atom for Z="<<z<<" is (cm^2)" 
187             <<sigma/cm/cm << G4endl;              179             <<sigma/cm/cm << G4endl;
188       G4cout<<"=== Cross section per atom for     180       G4cout<<"=== Cross section per atom for Z="<<z<<" is (cm^-1)=" 
189             <<sigma*atomicNDensity/(1./cm) <<     181             <<sigma*atomicNDensity/(1./cm) << G4endl;
190       G4cout<<"=== G4DNAQuinnPlasmonExcitation    182       G4cout<<"=== G4DNAQuinnPlasmonExcitationModel - XS INFO END" << G4endl;
191     }                                             183     }
192   }                                               184   } 
193                                                   185 
194   return sigma*atomicNDensity;                    186   return sigma*atomicNDensity;
195 }                                                 187 }
196                                                   188 
197 //....oooOO0OOooo........oooOO0OOooo........oo    189 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
198                                                   190 
199 void G4DNAQuinnPlasmonExcitationModel::SampleS    191 void G4DNAQuinnPlasmonExcitationModel::SampleSecondaries
200                          (std::vector<G4Dynami    192                          (std::vector<G4DynamicParticle*>* /*fvect*/,
201                           const G4MaterialCuts    193                           const G4MaterialCutsCouple* couple,
202                           const G4DynamicParti    194                           const G4DynamicParticle* aDynamicParticle,
203                           G4double,G4double)      195                           G4double,G4double)
204 {                                                 196 {
205                                                   197 
206   if (verboseLevel > 3)                           198   if (verboseLevel > 3)
207   {                                               199   {
208     G4cout <<                                     200     G4cout << 
209            "Calling SampleSecondaries() of G4D    201            "Calling SampleSecondaries() of G4DNAQuinnPlasmonExcitationModel"
210            << G4endl;                             202            << G4endl;
211   }                                               203   }
212                                                   204 
213   const G4Material *material = couple->GetMate    205   const G4Material *material = couple->GetMaterial();
214                                                   206 
215   G4ParticleDefinition* particle = aDynamicPar    207   G4ParticleDefinition* particle = aDynamicParticle->GetDefinition();
216                                                   208   
217   G4double k                     = aDynamicPar    209   G4double k                     = aDynamicParticle->GetKineticEnergy();
218                                                   210 
219   if(particle == G4Electron::ElectronDefinitio    211   if(particle == G4Electron::ElectronDefinition())
220   {                                               212   {
221     G4double  e      = 1.;                        213     G4double  e      = 1.;
222     G4int     z      = material->GetZ();          214     G4int     z      = material->GetZ();
223     G4int     Nve    = 0;                         215     G4int     Nve    = 0;
224                                                   216 
225     //TODO: have to be change to realistic!!      217     //TODO: have to be change to realistic!!
226     if(z<100) Nve    = nValenceElectron[z];       218     if(z<100) Nve    = nValenceElectron[z];
227                                                   219 
228     G4double  A      = material->GetA()/g/mole    220     G4double  A      = material->GetA()/g/mole;
229     G4double  Dens   = material->GetDensity()/    221     G4double  Dens   = material->GetDensity()/g*cm*cm*cm;
230     G4double  veDens = Dens*CLHEP::Avogadro*Nv    222     G4double  veDens = Dens*CLHEP::Avogadro*Nve/A;
231                                                   223 
232     G4double omega_p = std::sqrt(veDens*std::p    224     G4double omega_p = std::sqrt(veDens*std::pow(e,2)/
233                       (CLHEP::epsilon0/(1./cm)    225                       (CLHEP::epsilon0/(1./cm)*CLHEP::electron_mass_c2
234                       /(CLHEP::c_squared/cm/cm    226                       /(CLHEP::c_squared/cm/cm)));
235                                                   227     
236     G4double excitationEnergy   = CLHEP::hbar_    228     G4double excitationEnergy   = CLHEP::hbar_Planck*omega_p;
237     G4double newEnergy          = k - excitati    229     G4double newEnergy          = k - excitationEnergy;
238                                                   230 
239                                                   231 
240     if (newEnergy > 0)                            232     if (newEnergy > 0)
241     {                                             233     {
242       fParticleChangeForGamma->                   234       fParticleChangeForGamma->
243             ProposeMomentumDirection(aDynamicP    235             ProposeMomentumDirection(aDynamicParticle->GetMomentumDirection());
244                                                   236       
245       fParticleChangeForGamma->ProposeLocalEne    237       fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
246                                                   238       
247       if(!statCode)                               239       if(!statCode)
248       {                                           240       {
249            fParticleChangeForGamma->SetPropose    241            fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
250       }                                           242       }
251       else                                        243       else
252       {                                           244       {
253            fParticleChangeForGamma->SetPropose    245            fParticleChangeForGamma->SetProposedKineticEnergy(k);
254                                                   246       
255       }                                           247       }
256     }                                             248     }
257   }                                               249   } 
258 }                                                 250 }
259                                                   251 
260 //....oooOO0OOooo........oooOO0OOooo........oo    252 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
261                                                   253 
262 G4double G4DNAQuinnPlasmonExcitationModel::Get    254 G4double G4DNAQuinnPlasmonExcitationModel::GetCrossSection
263                          (const G4Material* ma    255                          (const G4Material* material,
264                           const G4ParticleDefi    256                           const G4ParticleDefinition* particle,
265                           G4double kineticEner    257                           G4double kineticEnergy)
266 {                                                 258 {
267   G4double value=0;                               259   G4double value=0; 
268                                                   260 
269   if(particle == G4Electron::ElectronDefinitio    261   if(particle == G4Electron::ElectronDefinition())
270   {                                               262   {
271      G4double  e      = 1.;                       263      G4double  e      = 1.;
272      G4int     z      = material->GetZ();         264      G4int     z      = material->GetZ();
273      G4int     Nve    = 0;                        265      G4int     Nve    = 0;
274      if(z<100) Nve    = nValenceElectron[z];      266      if(z<100) Nve    = nValenceElectron[z];
275      G4double  A      = material->GetA()/g/mol    267      G4double  A      = material->GetA()/g/mole;
276      G4double  Dens   = material->GetDensity()    268      G4double  Dens   = material->GetDensity()/g*cm*cm*cm;
277      G4double  veDens = Dens*CLHEP::Avogadro*N    269      G4double  veDens = Dens*CLHEP::Avogadro*Nve/A;
278                                                   270 
279      G4double omega_p = std::sqrt(veDens*std::    271      G4double omega_p = std::sqrt(veDens*std::pow(e,2)
280                         /(CLHEP::epsilon0/(1./    272                         /(CLHEP::epsilon0/(1./cm)*CLHEP::electron_mass_c2/
281                         (CLHEP::c_squared/cm/c    273                         (CLHEP::c_squared/cm/cm)));
282                                                   274      
283      G4double fEnergy = std::pow(CLHEP::h_Plan    275      G4double fEnergy = std::pow(CLHEP::h_Planck,2)/(8*CLHEP::electron_mass_c2)*
284                         std::pow(3*veDens/CLHE    276                         std::pow(3*veDens/CLHEP::pi,2./3.)/e
285                         *(CLHEP::c_squared/cm/    277                         *(CLHEP::c_squared/cm/cm);
286                                                   278      
287      G4double p0      = sqrt(2*CLHEP::electron    279      G4double p0      = sqrt(2*CLHEP::electron_mass_c2
288                         /(CLHEP::c_squared/cm/    280                         /(CLHEP::c_squared/cm/cm)*fEnergy);
289                                                   281      
290      G4double p       = sqrt(2*CLHEP::electron    282      G4double p       = sqrt(2*CLHEP::electron_mass_c2
291                         /(CLHEP::c_squared/cm/    283                         /(CLHEP::c_squared/cm/cm)*kineticEnergy);
292                                                   284 
293      G4double mfp     = 2*CLHEP::Bohr_radius/c    285      G4double mfp     = 2*CLHEP::Bohr_radius/cm*kineticEnergy
294                         /(CLHEP::hbar_Planck*o    286                         /(CLHEP::hbar_Planck*omega_p)/
295                         (G4Log((std::pow(std::    287                         (G4Log((std::pow(std::pow(p0,2)
296                         +2*CLHEP::electron_mas    288                         +2*CLHEP::electron_mass_c2/
297                         (CLHEP::c_squared/cm/c    289                         (CLHEP::c_squared/cm/cm)*omega_p
298                         *CLHEP::hbar_Planck,1.    290                         *CLHEP::hbar_Planck,1./2.)-p0)
299                         /(p-std::pow(std::pow(    291                         /(p-std::pow(std::pow(p,2)-2*CLHEP::electron_mass_c2/
300                         (CLHEP::c_squared/cm/c    292                         (CLHEP::c_squared/cm/cm)*omega_p
301                         *CLHEP::hbar_Planck,1.    293                         *CLHEP::hbar_Planck,1./2.))));
302                                                   294 
303      G4double excitationEnergy   = CLHEP::hbar    295      G4double excitationEnergy   = CLHEP::hbar_Planck*omega_p;
304                                                   296      
305      if((0<mfp)&&(0<veDens)&&(excitationEnergy    297      if((0<mfp)&&(0<veDens)&&(excitationEnergy<kineticEnergy)){
306        value            = 1./(veDens*mfp);        298        value            = 1./(veDens*mfp);
307      }                                            299      }
308   }                                               300   }
309   return value*cm*cm;                             301   return value*cm*cm;
310 }                                                 302 }
311                                                   303 
312 //....oooOO0OOooo........oooOO0OOooo........oo    304 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
313                                                   305 
314 G4int G4DNAQuinnPlasmonExcitationModel::GetNVa    306 G4int G4DNAQuinnPlasmonExcitationModel::GetNValenceElectron(G4int z)
315 {                                                 307 {
316                                                   308 
317   G4int Nve=0;                                    309   G4int Nve=0;
318                                                   310   
319   // Current limitation to gold                   311   // Current limitation to gold
320   if (z!=79){return 0.;}                          312   if (z!=79){return 0.;}
321                                                   313 
322   if (verboseLevel > 3)                           314   if (verboseLevel > 3)
323   {                                               315   {
324     G4cout <<                                     316     G4cout << 
325            "Calling GetNValenceElectron() of G    317            "Calling GetNValenceElectron() of G4DNAQuinnPlasmonExcitationModel"
326            << G4endl;                             318            << G4endl;
327   }                                               319   }
328                                                   320  
329   const char *datadir=nullptr;                 << 321   const char *datadir=0;
330                                                   322   
331   if(datadir == nullptr)                       << 323   if(!datadir)
332   {                                               324   {
333      datadir = G4FindDataDir("G4LEDATA");      << 325      datadir = getenv("G4LEDATA");
334      if(datadir == nullptr)                    << 326      if(!datadir)
335      {                                            327      {
336        G4Exception("G4DNAQuinnPlasmonExcitatio    328        G4Exception("G4DNAQuinnPlasmonExcitationModel::GetNValenceElectron()"
337                    ,"em0002",FatalException,      329                    ,"em0002",FatalException,
338                    "Enviroment variable G4LEDA    330                    "Enviroment variable G4LEDATA not defined");
339        return 0;                                  331        return 0;
340      }                                            332      }
341   }                                               333   }
342                                                   334 
343   std::ostringstream targetfile;                  335   std::ostringstream targetfile;
344   targetfile.str("");                             336   targetfile.str("");
345   targetfile.clear(stringstream::goodbit);        337   targetfile.clear(stringstream::goodbit);
346   targetfile << datadir <<"/dna/atomicstate_Z"    338   targetfile << datadir <<"/dna/atomicstate_Z"<< z <<".dat";
347   std::ifstream fin(targetfile.str().c_str());    339   std::ifstream fin(targetfile.str().c_str());
348                                                   340 
349   if(!fin)                                        341   if(!fin)
350   {                                               342   {
351     G4cout<< " Error : "<< targetfile.str() <<    343     G4cout<< " Error : "<< targetfile.str() <<" is not found "<<endl;
352     G4Exception("G4DNAQuinnPlasmonExcitationMo    344     G4Exception("G4DNAQuinnPlasmonExcitationModel::GetNValenceElectron()"
353                 ,"em0003",FatalException,         345                 ,"em0003",FatalException,
354                 "There is no target file");       346                 "There is no target file");
355     return 0;                                     347     return 0;
356   }                                               348   }
357                                                   349   
358   string buff0,buff1,buff2,buff3,buff4,buff5,b    350   string buff0,buff1,buff2,buff3,buff4,buff5,buff6;
359   fin >> buff0 >>buff1>>buff2>>buff3>>buff4>>b    351   fin >> buff0 >>buff1>>buff2>>buff3>>buff4>>buff5>>buff6;
360                                                   352   
361   while(true){                                    353   while(true){
362     fin >> buff0 >>buff1>>buff2>>buff3>>buff4>    354     fin >> buff0 >>buff1>>buff2>>buff3>>buff4>>buff5>>buff6;
363     if(!fin.eof())                                355     if(!fin.eof())
364     {                                             356     {
365       Nve = stoi(buff3);                          357       Nve = stoi(buff3);
366     }                                             358     }
367     else                                          359     else
368     {                                             360     {
369       break;                                      361       break;
370     }                                             362     }
371   }                                               363   }
372   return Nve;                                     364   return Nve;
373 }                                                 365 }
374                                                   366 
375                                                   367