Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/dna/models/src/G4DNASancheExcitationModel.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/G4DNASancheExcitationModel.cc (Version 11.3.0) and /processes/electromagnetic/dna/models/src/G4DNASancheExcitationModel.cc (Version 9.6.p2)


  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 // $Id$
 26 //                                                 27 //
 27                                                    28 
 28 // Created by Z. Francis                           29 // Created by Z. Francis
 29                                                    30 
 30 #include "G4DNASancheExcitationModel.hh"           31 #include "G4DNASancheExcitationModel.hh"
 31 #include "G4SystemOfUnits.hh"                      32 #include "G4SystemOfUnits.hh"
 32 #include "G4DNAMolecularMaterial.hh"               33 #include "G4DNAMolecularMaterial.hh"
 33                                                    34 
 34 //....oooOO0OOooo........oooOO0OOooo........oo <<  35 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 35                                                    36 
 36 using namespace std;                               37 using namespace std;
 37                                                    38 
 38 //#define SANCHE_VERBOSE                       <<  39 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 39                                                << 
 40 //....oooOO0OOooo........oooOO0OOooo........oo << 
 41                                                    40 
 42 G4DNASancheExcitationModel::G4DNASancheExcitat     41 G4DNASancheExcitationModel::G4DNASancheExcitationModel(const G4ParticleDefinition*,
 43                                                <<  42                                                        const G4String& nam)
 44     G4VEmModel(nam)                            <<  43     :G4VEmModel(nam),isInitialised(false)
 45 {                                                  44 {
 46   fpWaterDensity = nullptr;                    <<  45     //  nistwater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
                                                   >>  46     fpWaterDensity = 0;
 47                                                    47 
 48   SetLowEnergyLimit(2.*eV);                    <<  48     lowEnergyLimit = 2 * eV;
 49   SetHighEnergyLimit(100*eV);                  <<  49     highEnergyLimit = 100 * eV;
 50   nLevels = 9;                                 <<  50     SetLowEnergyLimit(lowEnergyLimit);
 51                                                <<  51     SetHighEnergyLimit(highEnergyLimit);
 52   verboseLevel = 0;                            <<  52     nLevels = 9;
 53   // Verbosity scale:                          <<  53 
 54   // 0 = nothing                               <<  54     verboseLevel= 0;
 55   // 1 = warning for energy non-conservation   <<  55     // Verbosity scale:
 56   // 2 = details of energy budget              <<  56     // 0 = nothing
 57   // 3 = calculation of cross sections, file o <<  57     // 1 = warning for energy non-conservation
 58   // 4 = entering in methods                   <<  58     // 2 = details of energy budget
                                                   >>  59     // 3 = calculation of cross sections, file openings, sampling of atoms
                                                   >>  60     // 4 = entering in methods
 59                                                    61 
 60 #ifdef SANCHE_VERBOSE                          <<  62     if (verboseLevel > 0)
 61   if (verboseLevel > 0)                        <<  63     {
 62   {                                            <<  64         G4cout << "Sanche Excitation model is constructed " << G4endl
 63     G4cout << "Sanche Excitation model is cons <<  65                << "Energy range: "
 64            << G4endl                           <<  66                << lowEnergyLimit / eV << " eV - "
 65            << "Energy range: "                 <<  67                << highEnergyLimit / eV << " eV"
 66            << LowEnergyLimit() / eV << " eV -  <<  68                << G4endl;
 67            << HighEnergyLimit() / eV << " eV"  <<  69     }
 68            << G4endl;                          <<  70     fParticleChangeForGamma = 0;
 69   }                                            <<  71     fpWaterDensity = 0;
 70 #endif                                         << 
 71                                                << 
 72   fParticleChangeForGamma = nullptr;           << 
 73   fpWaterDensity = nullptr;                    << 
 74                                                << 
 75   // Selection of stationary mode              << 
 76                                                << 
 77   statCode = false;                            << 
 78 }                                                  72 }
 79                                                    73 
 80 //....oooOO0OOooo........oooOO0OOooo........oo <<  74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 81                                                    75 
 82 G4DNASancheExcitationModel::~G4DNASancheExcita     76 G4DNASancheExcitationModel::~G4DNASancheExcitationModel()
 83 = default;                                     <<  77 {}
 84                                                    78 
 85 //....oooOO0OOooo........oooOO0OOooo........oo <<  79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 86                                                    80 
 87 void                                           <<  81 void G4DNASancheExcitationModel::Initialise(const G4ParticleDefinition* /*particle*/,
 88 G4DNASancheExcitationModel::                   <<  82                                             const G4DataVector& /*cuts*/)
 89 Initialise(const G4ParticleDefinition* /*parti << 
 90            const G4DataVector& /*cuts*/)       << 
 91 {                                                  83 {
 92                                                << 
 93 #ifdef SANCHE_VERBOSE                          << 
 94   if (verboseLevel > 3)                        << 
 95   {                                            << 
 96     G4cout << "Calling G4DNASancheExcitationMo << 
 97            << G4endl;                          << 
 98   }                                            << 
 99 #endif                                         << 
100                                                << 
101   // Energy limits                             << 
102                                                    84 
103   if (LowEnergyLimit() < 2.*eV)                << 
104   {                                            << 
105     G4Exception("*** WARNING : the G4DNASanche << 
106                 "validated below 2 eV !",      << 
107                 "", JustWarning, "");          << 
108   }                                            << 
109                                                    85 
110   if (HighEnergyLimit() > 100.*eV)             <<  86     if (verboseLevel > 3)
111   {                                            <<  87         G4cout << "Calling G4DNASancheExcitationModel::Initialise()" << G4endl;
112     G4cout << "G4DNASancheExcitationModel: hig << 
113     HighEnergyLimit()/eV << " eV to " << 100.  << 
114     SetHighEnergyLimit(100.*eV);               << 
115   }                                            << 
116                                                    88 
117   //                                           <<  89     // Energy limits
118 #ifdef SANCHE_VERBOSE                          << 
119   if (verboseLevel > 0)                        << 
120   {                                            << 
121     G4cout << "Sanche Excitation model is init << 
122     << "Energy range: "                        << 
123     << LowEnergyLimit() / eV << " eV - "       << 
124     << HighEnergyLimit() / eV << " eV"         << 
125     << G4endl;                                 << 
126   }                                            << 
127 #endif                                         << 
128                                                << 
129   // Initialize water density pointer          << 
130   fpWaterDensity = G4DNAMolecularMaterial::Ins << 
131       GetNumMolPerVolTableFor(G4Material::GetM << 
132                                                << 
133   if (isInitialised) {return;}                 << 
134                                                << 
135   fParticleChangeForGamma = GetParticleChangeF << 
136   isInitialised = true;                        << 
137                                                << 
138   const char *path = G4FindDataDir("G4LEDATA") << 
139   std::ostringstream eFullFileName;            << 
140   eFullFileName << path << "/dna/sigma_excitat << 
141   std::ifstream input(eFullFileName.str().c_st << 
142                                                    90 
143   if (!input)                                  <<  91     if (LowEnergyLimit() < lowEnergyLimit)
144   {                                            <<  92     {
145     G4Exception("G4DNASancheExcitationModel::I <<  93         G4cout << "G4DNASancheExcitationModel: low energy limit increased from " <<
146         FatalException,"Missing data file:/dna <<  94                   LowEnergyLimit()/eV << " eV to " << lowEnergyLimit/eV << " eV" << G4endl;
147   }                                            <<  95         SetLowEnergyLimit(lowEnergyLimit);
                                                   >>  96     }
148                                                    97 
149   // March 25th, 2014 - Vaclav Stepan, Sebasti <<  98     if (HighEnergyLimit() > highEnergyLimit)
150   // Added clear for MT                        <<  99     {
151   tdummyVec.clear();                           << 100         G4cout << "G4DNASancheExcitationModel: high energy limit decreased from " <<
152   //                                           << 101                   HighEnergyLimit()/eV << " eV to " << highEnergyLimit/eV << " eV" << G4endl;
                                                   >> 102         SetHighEnergyLimit(highEnergyLimit);
                                                   >> 103     }
153                                                   104 
154   G4double t;                                  << 105     //
155   G4double xs;                                 << 
156                                                   106 
157   while(!input.eof())                          << 107     if (verboseLevel > 0)
158   {                                            << 108     {
159     input>>t;                                  << 109         G4cout << "Sanche Excitation model is initialized " << G4endl
160     tdummyVec.push_back(t);                    << 110                << "Energy range: "
                                                   >> 111                << LowEnergyLimit() / eV << " eV - "
                                                   >> 112                << HighEnergyLimit() / eV << " eV"
                                                   >> 113                << G4endl;
                                                   >> 114     }
161                                                   115 
162     fEnergyLevelXS.emplace_back();             << 116     // Initialize water density pointer
163     fEnergyTotalXS.push_back(0);               << 117     fpWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
164     std::vector<G4double>& levelXS = fEnergyLe << 
165     levelXS.reserve(9);                        << 
166                                                   118 
167     // G4cout<<t;                              << 119     if (isInitialised) { return; }
                                                   >> 120     fParticleChangeForGamma = GetParticleChangeForGamma();
                                                   >> 121     isInitialised = true;
                                                   >> 122 
                                                   >> 123     char *path = getenv("G4LEDATA");
                                                   >> 124     std::ostringstream eFullFileName;
                                                   >> 125     eFullFileName << path << "/dna/sigma_excitationvib_e_sanche.dat";
                                                   >> 126     std::ifstream input(eFullFileName.str().c_str());
168                                                   127 
169     for(size_t i = 0 ; i < 9 ;++i)             << 128     if (!input)
170     {                                             129     {
171       input>>xs;                               << 130         G4Exception("G4DNASancheExcitationModel::Initialise","em0003",
172       levelXS.push_back(xs);                   << 131                     FatalException,"Missing data file:/dna/sigma_excitationvib_e_sanche.dat");
173       fEnergyTotalXS.back() += xs;             << 
174       // G4cout <<"  " << levelXS[i];          << 
175     }                                             132     }
176                                                   133 
177     // G4cout << G4endl;                       << 134     while(!input.eof())
178   }                                            << 135     {
                                                   >> 136         double t;
                                                   >> 137         input>>t;
                                                   >> 138         tdummyVec.push_back(t);
                                                   >> 139         input>>map1[t][0]>>map1[t][1]>>map1[t][2]>>map1[t][3]>>map1[t][4]>>map1[t][5]>>map1[t][6]>>map1[t][7]>>map1[t][8];
                                                   >> 140         //G4cout<<t<<"  "<<map1[t][0]<<map1[t][1]<<map1[t][2]<<map1[t][3]<<map1[t][4]<<map1[t][5]<<map1[t][6]<<map1[t][7]<<map1[t][8]<<G4endl;
                                                   >> 141     }
179 }                                                 142 }
180                                                   143 
181 //....oooOO0OOooo........oooOO0OOooo........oo << 144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
182                                                   145 
183 G4double G4DNASancheExcitationModel::CrossSect << 146 G4double G4DNASancheExcitationModel::CrossSectionPerVolume(const G4Material* material, 
184                                                << 147                                                            const G4ParticleDefinition* particleDefinition,
185 #ifdef SANCHE_VERBOSE                          << 
186                                                << 
187 #endif                                         << 
188                                                << 
189                                                   148                                                            G4double ekin,
190                                                   149                                                            G4double,
191                                                   150                                                            G4double)
192 {                                                 151 {
193 #ifdef SANCHE_VERBOSE                          << 152     if (verboseLevel > 3)
194   if (verboseLevel > 3)                        << 153         G4cout << "Calling CrossSectionPerVolume() of G4DNASancheExcitationModel" << G4endl;
195   {                                            << 
196     G4cout << "Calling CrossSectionPerVolume() << 
197            << G4endl;                          << 
198   }                                            << 
199 #endif                                         << 
200                                                   154 
201   // Calculate total cross section for model   << 155     // Calculate total cross section for model
202                                                   156 
203   G4double sigma = 0.;                         << 157     G4double sigma=0;
204                                                   158 
205   G4double waterDensity = (*fpWaterDensity)[ma << 159     G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
206                                                   160 
207   if (ekin >= LowEnergyLimit() && ekin <= High << 161     if(waterDensity!= 0.0)
208     sigma =  TotalCrossSection(ekin);          << 162         //  if (material == nistwater || material->GetBaseMaterial() == nistwater)
                                                   >> 163     {
209                                                   164 
210 #ifdef SANCHE_VERBOSE                          << 165         if (particleDefinition == G4Electron::ElectronDefinition())
211   if (verboseLevel > 2)                        << 166         {
212   {                                            << 167             if (ekin >= lowEnergyLimit && ekin < highEnergyLimit)
213     G4cout << "_______________________________ << 168             {
214     G4cout << "=== G4DNASancheExcitationModel  << 169                 sigma = Sum(ekin);
215     G4cout << "=== Kinetic energy(eV)=" << eki << 170             }
216     G4cout << "=== Cross section per water mol << 171         }
217     G4cout << "=== Cross section per water mol << 172 
218     G4cout << "=== G4DNASancheExcitationModel  << 173         if (verboseLevel > 2)
219   }                                            << 174         {
220 #endif                                         << 175             G4cout << "__________________________________" << G4endl;
                                                   >> 176             G4cout << "°°° G4DNASancheExcitationModel - XS INFO START" << G4endl;
                                                   >> 177             G4cout << "°°° Kinetic energy(eV)=" << ekin/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
                                                   >> 178             G4cout << "°°° Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
                                                   >> 179             G4cout << "°°° Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
                                                   >> 180             //      G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
                                                   >> 181             G4cout << "°°° G4DNASancheExcitationModel - XS INFO END" << G4endl;
                                                   >> 182         }
                                                   >> 183 
                                                   >> 184     } // if water
                                                   >> 185 
                                                   >> 186 
                                                   >> 187     //  return sigma*2*material->GetAtomicNumDensityVector()[1];
                                                   >> 188     return sigma*2*waterDensity;
                                                   >> 189     // see papers for factor 2 description
221                                                   190 
222   return sigma*2.*waterDensity;                << 
223   // see papers for factor 2 description (liqu << 
224 }                                                 191 }
225                                                   192 
226 //....oooOO0OOooo........oooOO0OOooo........oo << 193 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
227                                                   194 
228 void G4DNASancheExcitationModel::SampleSeconda << 195 void G4DNASancheExcitationModel::SampleSecondaries(std::vector<G4DynamicParticle*>*,
229                                                << 
230                                                   196                                                    const G4MaterialCutsCouple*,
231                                                   197                                                    const G4DynamicParticle* aDynamicElectron,
232                                                   198                                                    G4double,
233                                                   199                                                    G4double)
234 {                                                 200 {
235 #ifdef SANCHE_VERBOSE                          << 
236   if (verboseLevel > 3)                        << 
237   {                                            << 
238     G4cout << "Calling SampleSecondaries() of  << 
239            << G4endl;                          << 
240   }                                            << 
241 #endif                                         << 
242                                                   201 
243   G4double electronEnergy0 = aDynamicElectron- << 202     if (verboseLevel > 3)
244   G4int level = RandomSelect(electronEnergy0); << 203         G4cout << "Calling SampleSecondaries() of G4DNASancheExcitationModel" << G4endl;
245   G4double excitationEnergy = VibrationEnergy( << 
246   G4double newEnergy = electronEnergy0 - excit << 
247                                                << 
248   /*                                           << 
249    if (electronEnergy0 < highEnergyLimit)      << 
250    {                                           << 
251      if (newEnergy >= lowEnergyLimit)          << 
252      {                                         << 
253        fParticleChangeForGamma->ProposeMomentu << 
254        fParticleChangeForGamma->SetProposedKin << 
255        fParticleChangeForGamma->ProposeLocalEn << 
256      }                                         << 
257                                                << 
258      else                                      << 
259      {                                         << 
260        fParticleChangeForGamma->ProposeTrackSt << 
261        fParticleChangeForGamma->ProposeLocalEn << 
262      }                                         << 
263    }                                           << 
264    */                                          << 
265                                                   204 
266   if (electronEnergy0 <= HighEnergyLimit() &&  << 205     G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
267   {                                            << 206     G4int level = RandomSelect(electronEnergy0);
                                                   >> 207     G4double excitationEnergy = VibrationEnergy(level); // levels go from 0 to 8
                                                   >> 208     G4double newEnergy = electronEnergy0 - excitationEnergy;
268                                                   209 
269     if (!statCode)                             << 210     /*
270     {                                          << 211   if (electronEnergy0 < highEnergyLimit)
                                                   >> 212   {
                                                   >> 213     if (newEnergy >= lowEnergyLimit)
                                                   >> 214     {
271       fParticleChangeForGamma->ProposeMomentum    215       fParticleChangeForGamma->ProposeMomentumDirection(aDynamicElectron->GetMomentumDirection());
272       fParticleChangeForGamma->SetProposedKine    216       fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
273       fParticleChangeForGamma->ProposeLocalEne    217       fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
274     }                                             218     }
275                                                   219 
276     else                                       << 220     else
277     {                                             221     {
278       fParticleChangeForGamma->ProposeMomentum << 222       fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
279       fParticleChangeForGamma->SetProposedKine << 223       fParticleChangeForGamma->ProposeLocalEnergyDeposit(electronEnergy0);
280       fParticleChangeForGamma->ProposeLocalEne << 
281     }                                             224     }
282                                                << 
283   }                                               225   }
                                                   >> 226 */
284                                                   227 
285   //                                           << 228     if (electronEnergy0 < highEnergyLimit && newEnergy>0.)
286 }                                              << 229     {
287                                                << 230         fParticleChangeForGamma->ProposeMomentumDirection(aDynamicElectron->GetMomentumDirection());
288 //....oooOO0OOooo........oooOO0OOooo........oo << 231         fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
                                                   >> 232         fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
                                                   >> 233     }
289                                                   234 
290 G4double G4DNASancheExcitationModel::PartialCr << 235     //
291                                                << 
292 {                                              << 
293   // Protection against out of boundary access << 
294   if (t/eV==tdummyVec.back()) t=t*(1.-1e-12);  << 
295   //                                           << 
296                                                << 
297   auto t2 = std::upper_bound(tdummyVec.begin() << 
298                                                << 
299   auto t1 = t2 - 1;                            << 
300                                                << 
301   size_t i1 = t1 - tdummyVec.begin();          << 
302   size_t i2 = t2 - tdummyVec.begin();          << 
303                                                << 
304   G4double sigma = LinInterpolate((*t1), (*t2) << 
305                                 t / eV,        << 
306                                 fEnergyLevelXS << 
307                                 fEnergyLevelXS << 
308                                                << 
309   static const G4double conv_factor =  1e-16 * << 
310                                                << 
311   sigma *= conv_factor;                        << 
312   if (sigma == 0.) sigma = 1e-30;              << 
313   return (sigma);                              << 
314 }                                                 236 }
315                                                   237 
316 //....oooOO0OOooo........oooOO0OOooo........oo    238 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
317                                                   239 
318 G4double G4DNASancheExcitationModel::TotalCros << 240 G4double G4DNASancheExcitationModel::PartialCrossSection(G4double t, G4int level)
319 {                                                 241 {
320   // Protection against out of boundary access << 242     std::vector<double>::iterator t2 = std::upper_bound(tdummyVec.begin(),tdummyVec.end(), t/eV);
321   if (t/eV==tdummyVec.back()) t=t*(1.-1e-12);  << 243     std::vector<double>::iterator t1 = t2-1;
322   //                                           << 244 
323                                                << 245     double sigma = LinInterpolate((*t1), (*t2), t/eV, map1[*t1][level], map1[*t2][level]);
324   auto t2 = std::upper_bound(tdummyVec.begin() << 246     sigma*=1e-16*cm*cm;
325                                                << 247     if(sigma==0.)sigma=1e-30;
326   auto t1 = t2 - 1;                            << 248     return (sigma);
327                                                << 
328   size_t i1 = t1 - tdummyVec.begin();          << 
329   size_t i2 = t2 - tdummyVec.begin();          << 
330                                                << 
331   G4double sigma = LinInterpolate((*t1), (*t2) << 
332                                 t / eV,        << 
333                                 fEnergyTotalXS << 
334                                 fEnergyTotalXS << 
335                                                << 
336   static const G4double conv_factor =  1e-16 * << 
337                                                << 
338   sigma *= conv_factor;                        << 
339   if (sigma == 0.) sigma = 1e-30;              << 
340   return (sigma);                              << 
341 }                                                 249 }
342                                                   250 
343 //....oooOO0OOooo........oooOO0OOooo........oo    251 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
344                                                   252 
345 G4double G4DNASancheExcitationModel::Vibration    253 G4double G4DNASancheExcitationModel::VibrationEnergy(G4int level)
346 {                                                 254 {
347   static G4double energies[9] = { 0.01, 0.024, << 255     G4double energies[9] = {0.01, 0.024, 0.061, 0.092, 0.204, 0.417, 0.460, 0.500, 0.835};
348                            0.500, 0.835 };     << 256     return(energies[level]*eV);
349   return (energies[level] * eV);               << 
350 }                                                 257 }
351                                                   258 
352 //....oooOO0OOooo........oooOO0OOooo........oo    259 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
353                                                   260 
354 G4int G4DNASancheExcitationModel::RandomSelect    261 G4int G4DNASancheExcitationModel::RandomSelect(G4double k)
355 {                                                 262 {
356                                                   263 
357   // Level Selection Counting can be done here << 264     // Level Selection Counting can be done here !
358                                                   265 
359   G4int i = nLevels;                           << 266     G4int i = nLevels;
360   G4double value = 0.;                         << 267     G4double value = 0.;
361   std::deque<G4double> values;                 << 268     std::deque<double> values;
362                                                   269 
363   while (i > 0)                                << 270     while (i > 0)
364   {                                            << 271     {
365     i--;                                       << 272         i--;
366     G4double partial = PartialCrossSection(k,  << 273         G4double partial = PartialCrossSection(k,i);
367     values.push_front(partial);                << 274         values.push_front(partial);
368     value += partial;                          << 275         value += partial;
369   }                                            << 276     }
370                                                << 
371   value *= G4UniformRand();                    << 
372                                                   277 
373   i = nLevels;                                 << 278     value *= G4UniformRand();
                                                   >> 279     
                                                   >> 280     i = nLevels;
374                                                   281 
375   while (i > 0)                                << 282     while (i > 0)
376   {                                            << 
377     i--;                                       << 
378     if (values[i] > value)                     << 
379     {                                             283     {
380       //outcount<<i<<"  "<<VibrationEnergy(i)< << 284         i--;
381       return i;                                << 285         if (values[i] > value)
                                                   >> 286         {
                                                   >> 287             //outcount<<i<<"  "<<VibrationEnergy(i)<<G4endl;
                                                   >> 288             return i;
                                                   >> 289         }
                                                   >> 290         value -= values[i];
382     }                                             291     }
383     value -= values[i];                        << 
384   }                                            << 
385                                                   292 
386   //outcount<<0<<"  "<<VibrationEnergy(0)<<G4e << 293     //outcount<<0<<"  "<<VibrationEnergy(0)<<G4endl;
387                                                   294 
388   return 0;                                    << 295     return 0;
389 }                                                 296 }
390                                                   297 
391 //....oooOO0OOooo........oooOO0OOooo........oo    298 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
392                                                   299 
393 G4double G4DNASancheExcitationModel::Sum(G4dou    300 G4double G4DNASancheExcitationModel::Sum(G4double k)
394 {                                                 301 {
395   G4double totalCrossSection = 0.;             << 302     G4double totalCrossSection = 0.;
396                                                << 
397   for (G4int i = 0; i < nLevels; i++)          << 
398   {                                            << 
399     totalCrossSection += PartialCrossSection(k << 
400   }                                            << 
401                                                   303 
402   return totalCrossSection;                    << 304     for (G4int i=0; i<nLevels; i++)
                                                   >> 305     {
                                                   >> 306         totalCrossSection += PartialCrossSection(k,i);
                                                   >> 307     }
                                                   >> 308     return totalCrossSection;
403 }                                                 309 }
404                                                   310 
405 //....oooOO0OOooo........oooOO0OOooo........oo    311 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
406                                                   312 
407 G4double G4DNASancheExcitationModel::LinInterp << 313 G4double G4DNASancheExcitationModel::LinInterpolate(G4double e1, 
408                                                   314                                                     G4double e2,
409                                                   315                                                     G4double e,
410                                                   316                                                     G4double xs1,
411                                                   317                                                     G4double xs2)
412 {                                                 318 {
413   G4double a = (xs2 - xs1) / (e2 - e1);        << 319     G4double a = (xs2 - xs1) / (e2 - e1);
414   G4double b = xs2 - a * e2;                   << 320     G4double b = xs2 - a*e2;
415   G4double value = a * e + b;                  << 321     G4double value = a*e + b;
416   // G4cout<<"interP >>  "<<e1<<"  "<<e2<<"  " << 322     // G4cout<<"interP >>  "<<e1<<"  "<<e2<<"  "<<e<<"  "<<xs1<<"  "<<xs2<<"  "<<a<<"  "<<b<<"  "<<value<<G4endl;
417   // <<xs1<<"  "<<xs2<<"  "<<a<<"  "<<b<<"  "< << 
418                                                   323 
419   return value;                                << 324     return value;
420 }                                                 325 }
421                                                   326 
422                                                   327