G4DNASancheExcitationModel.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028 
00029 // Created by Z. Francis
00030 
00031 #include "G4DNASancheExcitationModel.hh"
00032 #include "G4SystemOfUnits.hh"
00033 #include "G4DNAMolecularMaterial.hh"
00034 
00035 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00036 
00037 using namespace std;
00038 
00039 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00040 
00041 G4DNASancheExcitationModel::G4DNASancheExcitationModel(const G4ParticleDefinition*,
00042                                                        const G4String& nam)
00043     :G4VEmModel(nam),isInitialised(false)
00044 {
00045     //  nistwater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
00046     fpWaterDensity = 0;
00047 
00048     lowEnergyLimit = 2 * eV;
00049     highEnergyLimit = 100 * eV;
00050     SetLowEnergyLimit(lowEnergyLimit);
00051     SetHighEnergyLimit(highEnergyLimit);
00052     nLevels = 9;
00053 
00054     verboseLevel= 0;
00055     // Verbosity scale:
00056     // 0 = nothing
00057     // 1 = warning for energy non-conservation
00058     // 2 = details of energy budget
00059     // 3 = calculation of cross sections, file openings, sampling of atoms
00060     // 4 = entering in methods
00061 
00062     if (verboseLevel > 0)
00063     {
00064         G4cout << "Sanche Excitation model is constructed " << G4endl
00065                << "Energy range: "
00066                << lowEnergyLimit / eV << " eV - "
00067                << highEnergyLimit / eV << " eV"
00068                << G4endl;
00069     }
00070     fParticleChangeForGamma = 0;
00071     fpWaterDensity = 0;
00072 }
00073 
00074 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00075 
00076 G4DNASancheExcitationModel::~G4DNASancheExcitationModel()
00077 {}
00078 
00079 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00080 
00081 void G4DNASancheExcitationModel::Initialise(const G4ParticleDefinition* /*particle*/,
00082                                             const G4DataVector& /*cuts*/)
00083 {
00084 
00085 
00086     if (verboseLevel > 3)
00087         G4cout << "Calling G4DNASancheExcitationModel::Initialise()" << G4endl;
00088 
00089     // Energy limits
00090 
00091     if (LowEnergyLimit() < lowEnergyLimit)
00092     {
00093         G4cout << "G4DNASancheExcitationModel: low energy limit increased from " <<
00094                   LowEnergyLimit()/eV << " eV to " << lowEnergyLimit/eV << " eV" << G4endl;
00095         SetLowEnergyLimit(lowEnergyLimit);
00096     }
00097 
00098     if (HighEnergyLimit() > highEnergyLimit)
00099     {
00100         G4cout << "G4DNASancheExcitationModel: high energy limit decreased from " <<
00101                   HighEnergyLimit()/eV << " eV to " << highEnergyLimit/eV << " eV" << G4endl;
00102         SetHighEnergyLimit(highEnergyLimit);
00103     }
00104 
00105     //
00106 
00107     if (verboseLevel > 0)
00108     {
00109         G4cout << "Sanche Excitation model is initialized " << G4endl
00110                << "Energy range: "
00111                << LowEnergyLimit() / eV << " eV - "
00112                << HighEnergyLimit() / eV << " eV"
00113                << G4endl;
00114     }
00115 
00116     // Initialize water density pointer
00117     fpWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
00118 
00119     if (isInitialised) { return; }
00120     fParticleChangeForGamma = GetParticleChangeForGamma();
00121     isInitialised = true;
00122 
00123     char *path = getenv("G4LEDATA");
00124     std::ostringstream eFullFileName;
00125     eFullFileName << path << "/dna/sigma_excitationvib_e_sanche.dat";
00126     std::ifstream input(eFullFileName.str().c_str());
00127 
00128     if (!input)
00129     {
00130         G4Exception("G4DNASancheExcitationModel::Initialise","em0003",
00131                     FatalException,"Missing data file:/dna/sigma_excitationvib_e_sanche.dat");
00132     }
00133 
00134     while(!input.eof())
00135     {
00136         double t;
00137         input>>t;
00138         tdummyVec.push_back(t);
00139         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];
00140         //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;
00141     }
00142 }
00143 
00144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00145 
00146 G4double G4DNASancheExcitationModel::CrossSectionPerVolume(const G4Material* material, 
00147                                                            const G4ParticleDefinition* particleDefinition,
00148                                                            G4double ekin,
00149                                                            G4double,
00150                                                            G4double)
00151 {
00152     if (verboseLevel > 3)
00153         G4cout << "Calling CrossSectionPerVolume() of G4DNASancheExcitationModel" << G4endl;
00154 
00155     // Calculate total cross section for model
00156 
00157     G4double sigma=0;
00158 
00159     G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
00160 
00161     if(waterDensity!= 0.0)
00162         //  if (material == nistwater || material->GetBaseMaterial() == nistwater)
00163     {
00164 
00165         if (particleDefinition == G4Electron::ElectronDefinition())
00166         {
00167             if (ekin >= lowEnergyLimit && ekin < highEnergyLimit)
00168             {
00169                 sigma = Sum(ekin);
00170             }
00171         }
00172 
00173         if (verboseLevel > 2)
00174         {
00175             G4cout << "__________________________________" << G4endl;
00176             G4cout << "°°° G4DNASancheExcitationModel - XS INFO START" << G4endl;
00177             G4cout << "°°° Kinetic energy(eV)=" << ekin/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
00178             G4cout << "°°° Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
00179             G4cout << "°°° Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
00180             //      G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
00181             G4cout << "°°° G4DNASancheExcitationModel - XS INFO END" << G4endl;
00182         }
00183 
00184     } // if water
00185 
00186 
00187     //  return sigma*2*material->GetAtomicNumDensityVector()[1];
00188     return sigma*2*waterDensity;
00189     // see papers for factor 2 description
00190 
00191 }
00192 
00193 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00194 
00195 void G4DNASancheExcitationModel::SampleSecondaries(std::vector<G4DynamicParticle*>*,
00196                                                    const G4MaterialCutsCouple*,
00197                                                    const G4DynamicParticle* aDynamicElectron,
00198                                                    G4double,
00199                                                    G4double)
00200 {
00201 
00202     if (verboseLevel > 3)
00203         G4cout << "Calling SampleSecondaries() of G4DNASancheExcitationModel" << G4endl;
00204 
00205     G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
00206     G4int level = RandomSelect(electronEnergy0);
00207     G4double excitationEnergy = VibrationEnergy(level); // levels go from 0 to 8
00208     G4double newEnergy = electronEnergy0 - excitationEnergy;
00209 
00210     /*
00211   if (electronEnergy0 < highEnergyLimit)
00212   {
00213     if (newEnergy >= lowEnergyLimit)
00214     {
00215       fParticleChangeForGamma->ProposeMomentumDirection(aDynamicElectron->GetMomentumDirection());
00216       fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
00217       fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
00218     }
00219 
00220     else
00221     {
00222       fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
00223       fParticleChangeForGamma->ProposeLocalEnergyDeposit(electronEnergy0);
00224     }
00225   }
00226 */
00227 
00228     if (electronEnergy0 < highEnergyLimit && newEnergy>0.)
00229     {
00230         fParticleChangeForGamma->ProposeMomentumDirection(aDynamicElectron->GetMomentumDirection());
00231         fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
00232         fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
00233     }
00234 
00235     //
00236 }
00237 
00238 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00239 
00240 G4double G4DNASancheExcitationModel::PartialCrossSection(G4double t, G4int level)
00241 {
00242     std::vector<double>::iterator t2 = std::upper_bound(tdummyVec.begin(),tdummyVec.end(), t/eV);
00243     std::vector<double>::iterator t1 = t2-1;
00244 
00245     double sigma = LinInterpolate((*t1), (*t2), t/eV, map1[*t1][level], map1[*t2][level]);
00246     sigma*=1e-16*cm*cm;
00247     if(sigma==0.)sigma=1e-30;
00248     return (sigma);
00249 }
00250 
00251 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00252 
00253 G4double G4DNASancheExcitationModel::VibrationEnergy(G4int level)
00254 {
00255     G4double energies[9] = {0.01, 0.024, 0.061, 0.092, 0.204, 0.417, 0.460, 0.500, 0.835};
00256     return(energies[level]*eV);
00257 }
00258 
00259 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00260 
00261 G4int G4DNASancheExcitationModel::RandomSelect(G4double k)
00262 {
00263 
00264     // Level Selection Counting can be done here !
00265 
00266     G4int i = nLevels;
00267     G4double value = 0.;
00268     std::deque<double> values;
00269 
00270     while (i > 0)
00271     {
00272         i--;
00273         G4double partial = PartialCrossSection(k,i);
00274         values.push_front(partial);
00275         value += partial;
00276     }
00277 
00278     value *= G4UniformRand();
00279     
00280     i = nLevels;
00281 
00282     while (i > 0)
00283     {
00284         i--;
00285         if (values[i] > value)
00286         {
00287             //outcount<<i<<"  "<<VibrationEnergy(i)<<G4endl;
00288             return i;
00289         }
00290         value -= values[i];
00291     }
00292 
00293     //outcount<<0<<"  "<<VibrationEnergy(0)<<G4endl;
00294 
00295     return 0;
00296 }
00297 
00298 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00299 
00300 G4double G4DNASancheExcitationModel::Sum(G4double k)
00301 {
00302     G4double totalCrossSection = 0.;
00303 
00304     for (G4int i=0; i<nLevels; i++)
00305     {
00306         totalCrossSection += PartialCrossSection(k,i);
00307     }
00308     return totalCrossSection;
00309 }
00310 
00311 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00312 
00313 G4double G4DNASancheExcitationModel::LinInterpolate(G4double e1, 
00314                                                     G4double e2,
00315                                                     G4double e,
00316                                                     G4double xs1,
00317                                                     G4double xs2)
00318 {
00319     G4double a = (xs2 - xs1) / (e2 - e1);
00320     G4double b = xs2 - a*e2;
00321     G4double value = a*e + b;
00322     // G4cout<<"interP >>  "<<e1<<"  "<<e2<<"  "<<e<<"  "<<xs1<<"  "<<xs2<<"  "<<a<<"  "<<b<<"  "<<value<<G4endl;
00323 
00324     return value;
00325 }
00326 

Generated on Mon May 27 17:48:03 2013 for Geant4 by  doxygen 1.4.7