G4DNARuddIonisationExtendedModel.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 // GEANT4 tag $Name:  $
00028 //
00029 
00030 
00031 // Modified by Z. Francis to handle HZE && inverse rudd function sampling 26-10-2010
00032 
00033 #include "G4DNARuddIonisationExtendedModel.hh"
00034 #include "G4PhysicalConstants.hh"
00035 #include "G4SystemOfUnits.hh"
00036 #include "G4UAtomicDeexcitation.hh"
00037 #include "G4LossTableManager.hh"
00038 #include "G4DNAChemistryManager.hh"
00039 #include "G4DNAMolecularMaterial.hh"
00040 
00041 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00042 
00043 using namespace std;
00044 
00045 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00046 
00047 G4DNARuddIonisationExtendedModel::G4DNARuddIonisationExtendedModel(const G4ParticleDefinition*,
00048                                                                    const G4String& nam)
00049     :G4VEmModel(nam),isInitialised(false)
00050 {
00051     //  nistwater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
00052     fpWaterDensity = 0;
00053 
00054     slaterEffectiveCharge[0]=0.;
00055     slaterEffectiveCharge[1]=0.;
00056     slaterEffectiveCharge[2]=0.;
00057     sCoefficient[0]=0.;
00058     sCoefficient[1]=0.;
00059     sCoefficient[2]=0.;
00060 
00061     lowEnergyLimitForA[1] = 0 * eV;
00062     lowEnergyLimitForA[2] = 0 * eV;
00063     lowEnergyLimitForA[3] = 0 * eV;
00064     lowEnergyLimitOfModelForA[1] = 100 * eV;
00065     lowEnergyLimitOfModelForA[4] = 1 * keV;
00066     lowEnergyLimitOfModelForA[5] = 0.5 * MeV; // For A = 3 or above, limit is MeV/uma
00067     killBelowEnergyForA[1] = lowEnergyLimitOfModelForA[1];
00068     killBelowEnergyForA[4] = lowEnergyLimitOfModelForA[4];
00069     killBelowEnergyForA[5] = lowEnergyLimitOfModelForA[5];
00070 
00071 
00072     verboseLevel= 0;
00073     // Verbosity scale:
00074     // 0 = nothing
00075     // 1 = warning for energy non-conservation
00076     // 2 = details of energy budget
00077     // 3 = calculation of cross sections, file openings, sampling of atoms
00078     // 4 = entering in methods
00079 
00080     if( verboseLevel>0 )
00081     {
00082         G4cout << "Rudd ionisation model is constructed " << G4endl;
00083     }
00084 
00085     //Mark this model as "applicable" for atomic deexcitation
00086     SetDeexcitationFlag(true);
00087     fAtomDeexcitation = 0;
00088     fParticleChangeForGamma = 0;
00089 }
00090 
00091 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00092 
00093 G4DNARuddIonisationExtendedModel::~G4DNARuddIonisationExtendedModel()
00094 {  
00095     // Cross section
00096 
00097     std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
00098     for (pos = tableData.begin(); pos != tableData.end(); ++pos)
00099     {
00100         G4DNACrossSectionDataSet* table = pos->second;
00101         delete table;
00102     }
00103 
00104     // The following removal is forbidden G4VEnergyLossModel takes care of deletion
00105     // however coverity will signal this as an error
00106     //if (fAtomDeexcitation) {delete  fAtomDeexcitation;}
00107 
00108 }
00109 
00110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00111 
00112 void G4DNARuddIonisationExtendedModel::Initialise(const G4ParticleDefinition* particle,
00113                                                   const G4DataVector& /*cuts*/)
00114 {
00115 
00116     if (verboseLevel > 3)
00117         G4cout << "Calling G4DNARuddIonisationExtendedModel::Initialise()" << G4endl;
00118 
00119     // Energy limits
00120 
00121     G4String fileProton("dna/sigma_ionisation_p_rudd");
00122     G4String fileHydrogen("dna/sigma_ionisation_h_rudd");
00123     G4String fileAlphaPlusPlus("dna/sigma_ionisation_alphaplusplus_rudd");
00124     G4String fileAlphaPlus("dna/sigma_ionisation_alphaplus_rudd");
00125     G4String fileHelium("dna/sigma_ionisation_he_rudd");
00126     G4String fileCarbon("dna/sigma_ionisation_c_rudd");
00127     G4String fileNitrogen("dna/sigma_ionisation_n_rudd");
00128     G4String fileOxygen("dna/sigma_ionisation_o_rudd");
00129     G4String fileIron("dna/sigma_ionisation_fe_rudd");
00130 
00131     G4DNAGenericIonsManager *instance;
00132     instance = G4DNAGenericIonsManager::Instance();
00133     G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
00134     G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
00135     G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
00136     G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
00137     G4ParticleDefinition* heliumDef = instance->GetIon("helium");
00138     G4ParticleDefinition* carbonDef = instance->GetIon("carbon");
00139     G4ParticleDefinition* nitrogenDef = instance->GetIon("nitrogen");
00140     G4ParticleDefinition* oxygenDef = instance->GetIon("oxygen");
00141     G4ParticleDefinition* ironDef = instance->GetIon("iron");
00142 
00143     G4String proton;
00144     G4String hydrogen;
00145     G4String alphaPlusPlus;
00146     G4String alphaPlus;
00147     G4String helium;
00148     G4String carbon;
00149     G4String nitrogen;
00150     G4String oxygen;
00151     G4String iron;
00152 
00153     G4double scaleFactor = 1 * m*m;
00154 
00155     // LIMITS AND DATA
00156 
00157     proton = protonDef->GetParticleName();
00158     tableFile[proton] = fileProton;
00159     lowEnergyLimit[proton] = lowEnergyLimitForA[1];
00160     highEnergyLimit[proton] = 500. * keV;
00161 
00162     // Cross section
00163 
00164     G4DNACrossSectionDataSet* tableProton = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
00165                                                                          eV,
00166                                                                          scaleFactor );
00167     tableProton->LoadData(fileProton);
00168     tableData[proton] = tableProton;
00169 
00170     // **********************************************************************************************
00171     
00172     hydrogen = hydrogenDef->GetParticleName();
00173     tableFile[hydrogen] = fileHydrogen;
00174 
00175     lowEnergyLimit[hydrogen] = lowEnergyLimitForA[1];
00176     highEnergyLimit[hydrogen] = 100. * MeV;
00177 
00178     // Cross section
00179 
00180     G4DNACrossSectionDataSet* tableHydrogen = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
00181                                                                            eV,
00182                                                                            scaleFactor );
00183     tableHydrogen->LoadData(fileHydrogen);
00184 
00185     tableData[hydrogen] = tableHydrogen;
00186 
00187     // **********************************************************************************************
00188     
00189     alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
00190     tableFile[alphaPlusPlus] = fileAlphaPlusPlus;
00191 
00192     lowEnergyLimit[alphaPlusPlus] = lowEnergyLimitForA[4];
00193     highEnergyLimit[alphaPlusPlus] = 400. * MeV;
00194 
00195     // Cross section
00196 
00197     G4DNACrossSectionDataSet* tableAlphaPlusPlus = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
00198                                                                                 eV,
00199                                                                                 scaleFactor );
00200     tableAlphaPlusPlus->LoadData(fileAlphaPlusPlus);
00201 
00202     tableData[alphaPlusPlus] = tableAlphaPlusPlus;
00203 
00204     // **********************************************************************************************
00205     
00206     alphaPlus = alphaPlusDef->GetParticleName();
00207     tableFile[alphaPlus] = fileAlphaPlus;
00208 
00209     lowEnergyLimit[alphaPlus] = lowEnergyLimitForA[4];
00210     highEnergyLimit[alphaPlus] = 400. * MeV;
00211 
00212     // Cross section
00213 
00214     G4DNACrossSectionDataSet* tableAlphaPlus = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
00215                                                                             eV,
00216                                                                             scaleFactor );
00217     tableAlphaPlus->LoadData(fileAlphaPlus);
00218     tableData[alphaPlus] = tableAlphaPlus;
00219 
00220     // **********************************************************************************************
00221     
00222     helium = heliumDef->GetParticleName();
00223     tableFile[helium] = fileHelium;
00224 
00225     lowEnergyLimit[helium] = lowEnergyLimitForA[4];
00226     highEnergyLimit[helium] = 400. * MeV;
00227 
00228     // Cross section
00229 
00230     G4DNACrossSectionDataSet* tableHelium = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
00231                                                                          eV,
00232                                                                          scaleFactor );
00233     tableHelium->LoadData(fileHelium);
00234     tableData[helium] = tableHelium;
00235 
00236     // **********************************************************************************************
00237     
00238     carbon = carbonDef->GetParticleName();
00239     tableFile[carbon] = fileCarbon;
00240 
00241     lowEnergyLimit[carbon] = lowEnergyLimitForA[5] * particle->GetAtomicMass();
00242     highEnergyLimit[carbon] = 1e6* particle->GetAtomicMass() * MeV;
00243 
00244     // Cross section
00245 
00246     G4DNACrossSectionDataSet* tableCarbon = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
00247                                                                          eV,
00248                                                                          scaleFactor );
00249     tableCarbon->LoadData(fileCarbon);
00250     tableData[carbon] = tableCarbon;
00251 
00252     // **********************************************************************************************
00253     
00254     oxygen = oxygenDef->GetParticleName();
00255     tableFile[oxygen] = fileOxygen;
00256 
00257     lowEnergyLimit[oxygen] = lowEnergyLimitForA[5]* particle->GetAtomicMass();
00258     highEnergyLimit[oxygen] = 1e6* particle->GetAtomicMass()* MeV;
00259 
00260     // Cross section
00261 
00262     G4DNACrossSectionDataSet* tableOxygen = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
00263                                                                          eV,
00264                                                                          scaleFactor );
00265     tableOxygen->LoadData(fileOxygen);
00266     tableData[oxygen] = tableOxygen;
00267 
00268     // **********************************************************************************************
00269     
00270     nitrogen = nitrogenDef->GetParticleName();
00271     tableFile[nitrogen] = fileNitrogen;
00272 
00273     lowEnergyLimit[nitrogen] = lowEnergyLimitForA[5]* particle->GetAtomicMass();
00274     highEnergyLimit[nitrogen] = 1e6* particle->GetAtomicMass()* MeV;
00275 
00276     // Cross section
00277 
00278     G4DNACrossSectionDataSet* tableNitrogen = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
00279                                                                            eV,
00280                                                                            scaleFactor );
00281     tableNitrogen->LoadData(fileNitrogen);
00282     tableData[nitrogen] = tableNitrogen;
00283 
00284     // **********************************************************************************************
00285     
00286     iron = ironDef->GetParticleName();
00287     tableFile[iron] = fileIron;
00288 
00289     lowEnergyLimit[iron] = lowEnergyLimitForA[5]* particle->GetAtomicMass();
00290     highEnergyLimit[iron] = 1e6* particle->GetAtomicMass()* MeV;
00291 
00292     // Cross section
00293 
00294     G4DNACrossSectionDataSet* tableIron = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
00295                                                                        eV,
00296                                                                        scaleFactor );
00297     tableIron->LoadData(fileIron);
00298     tableData[iron] = tableIron;
00299 
00300     // ZF Following lines can be replaced by:
00301     SetLowEnergyLimit(lowEnergyLimit[particle->GetParticleName()]);
00302     SetHighEnergyLimit(highEnergyLimit[particle->GetParticleName()]);
00303     // at least for HZE
00304 
00305     /*
00306   if (particle==protonDef)
00307   {
00308     SetLowEnergyLimit(lowEnergyLimit[proton]);
00309     SetHighEnergyLimit(highEnergyLimit[proton]);
00310   }
00311 
00312   if (particle==hydrogenDef)
00313   {
00314     SetLowEnergyLimit(lowEnergyLimit[hydrogen]);
00315     SetHighEnergyLimit(highEnergyLimit[hydrogen]);
00316   }
00317 
00318   if (particle==heliumDef)
00319   {
00320     SetLowEnergyLimit(lowEnergyLimit[helium]);
00321     SetHighEnergyLimit(highEnergyLimit[helium]);
00322   }
00323 
00324   if (particle==alphaPlusDef)
00325   {
00326     SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
00327     SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
00328   }
00329 
00330   if (particle==alphaPlusPlusDef)
00331   {
00332     SetLowEnergyLimit(lowEnergyLimit[alphaPlusPlus]);
00333     SetHighEnergyLimit(highEnergyLimit[alphaPlusPlus]);
00334   }
00335 
00336   if (particle==carbonDef)
00337   {
00338     SetLowEnergyLimit(lowEnergyLimit[carbon]);
00339     SetHighEnergyLimit(highEnergyLimit[carbon]);
00340   }
00341 
00342   if (particle==oxygenDef)
00343   {
00344     SetLowEnergyLimit(lowEnergyLimit[oxygen]);
00345     SetHighEnergyLimit(highEnergyLimit[oxygen]);
00346   }*/
00347 
00348     //----------------------------------------------------------------------
00349 
00350     if( verboseLevel>0 )
00351     {
00352         G4cout << "Rudd ionisation model is initialized " << G4endl
00353                << "Energy range: "
00354                << LowEnergyLimit() / eV << " eV - "
00355                << HighEnergyLimit() / keV << " keV for "
00356                << particle->GetParticleName()
00357                << G4endl;
00358     }
00359 
00360     // Initialize water density pointer
00361     fpWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
00362 
00363     //
00364 
00365     fAtomDeexcitation  = G4LossTableManager::Instance()->AtomDeexcitation();
00366 
00367     if (isInitialised) { return; }
00368     fParticleChangeForGamma = GetParticleChangeForGamma();
00369     isInitialised = true;
00370 }
00371 
00372 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00373 
00374 G4double G4DNARuddIonisationExtendedModel::CrossSectionPerVolume(const G4Material* material,
00375                                                                  const G4ParticleDefinition* particleDefinition,
00376                                                                  G4double k,
00377                                                                  G4double,
00378                                                                  G4double)
00379 {
00380     if (verboseLevel > 3)
00381         G4cout << "Calling CrossSectionPerVolume() of G4DNARuddIonisationExtendedModel" << G4endl;
00382 
00383     // Calculate total cross section for model
00384 
00385     G4DNAGenericIonsManager *instance;
00386     instance = G4DNAGenericIonsManager::Instance();
00387 
00388     if (
00389             particleDefinition != G4Proton::ProtonDefinition()
00390             &&
00391             particleDefinition != instance->GetIon("hydrogen")
00392             &&
00393             particleDefinition != instance->GetIon("alpha++")
00394             &&
00395             particleDefinition != instance->GetIon("alpha+")
00396             &&
00397             particleDefinition != instance->GetIon("helium")
00398             &&
00399             particleDefinition != instance->GetIon("carbon")
00400             &&
00401             particleDefinition != instance->GetIon("nitrogen")
00402             &&
00403             particleDefinition != instance->GetIon("oxygen")
00404             &&
00405             particleDefinition != instance->GetIon("iron")
00406             )
00407 
00408         return 0;
00409 
00410     G4double lowLim = 0;
00411 
00412     if (     particleDefinition == G4Proton::ProtonDefinition()
00413              ||  particleDefinition == instance->GetIon("hydrogen")
00414              )
00415 
00416         lowLim = lowEnergyLimitOfModelForA[1];
00417 
00418     else if (     particleDefinition == instance->GetIon("alpha++")
00419                   ||       particleDefinition == instance->GetIon("alpha+")
00420                   ||       particleDefinition == instance->GetIon("helium")
00421                   )
00422 
00423         lowLim = lowEnergyLimitOfModelForA[4];
00424 
00425     else lowLim = lowEnergyLimitOfModelForA[5];
00426 
00427     G4double highLim = 0;
00428     G4double sigma=0;
00429 
00430 
00431     G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
00432 
00433     if(waterDensity!= 0.0)
00434 //    if (material == nistwater || material->GetBaseMaterial() == nistwater)
00435     {
00436         const G4String& particleName = particleDefinition->GetParticleName();
00437 
00438         std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
00439         pos2 = highEnergyLimit.find(particleName);
00440 
00441         if (pos2 != highEnergyLimit.end())
00442         {
00443             highLim = pos2->second;
00444         }
00445 
00446         if (k <= highLim)
00447         {
00448 
00449             //SI : XS must not be zero otherwise sampling of secondaries method ignored
00450 
00451             if (k < lowLim) k = lowLim;
00452 
00453             //
00454 
00455             std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
00456             pos = tableData.find(particleName);
00457 
00458             if (pos != tableData.end())
00459             {
00460                 G4DNACrossSectionDataSet* table = pos->second;
00461                 if (table != 0)
00462                 {
00463                     sigma = table->FindValue(k);
00464                 }
00465             }
00466             else
00467             {
00468                 G4Exception("G4DNARuddIonisationExtendedModel::CrossSectionPerVolume","em0002",
00469                             FatalException,"Model not applicable to particle type.");
00470             }
00471 
00472         } // if (k >= lowLim && k < highLim)
00473 
00474         if (verboseLevel > 2)
00475         {
00476             G4cout << "__________________________________" << G4endl;
00477             G4cout << "°°° G4DNARuddIonisationExtendedModel - XS INFO START" << G4endl;
00478             G4cout << "°°° Kinetic energy(eV)=" << k/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
00479             G4cout << "°°° Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
00480             G4cout << "°°° Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
00481             //      G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
00482             G4cout << "°°° G4DNARuddIonisationExtendedModel - XS INFO END" << G4endl;
00483 
00484         }
00485 
00486     } // if (waterMaterial)
00487 
00488     return sigma*waterDensity;
00489 //    return sigma*material->GetAtomicNumDensityVector()[1];
00490 
00491 }
00492 
00493 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00494 
00495 void G4DNARuddIonisationExtendedModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
00496                                                          const G4MaterialCutsCouple* /*couple*/,
00497                                                          const G4DynamicParticle* particle,
00498                                                          G4double,
00499                                                          G4double)
00500 {
00501     if (verboseLevel > 3)
00502         G4cout << "Calling SampleSecondaries() of G4DNARuddIonisationExtendedModel" << G4endl;
00503 
00504     G4double lowLim = 0;
00505     G4double highLim = 0;
00506 
00507     // ZF: the following line summarizes the commented part
00508     if(particle->GetDefinition()->GetAtomicMass() <= 4) lowLim = killBelowEnergyForA[particle->GetDefinition()->GetAtomicMass()];
00509     else lowLim = killBelowEnergyForA[5]*particle->GetDefinition()->GetAtomicMass();
00510 
00511     /*
00512   if(particle->GetDefinition()->GetAtomicMass() >= 5) lowLim = killBelowEnergyForA[5]*particle->GetDefinition()->GetAtomicMass();
00513 
00514 
00515   if (     particle->GetDefinition() == G4Proton::ProtonDefinition()
00516        ||  particle->GetDefinition() == instance->GetIon("hydrogen")
00517      )
00518 
00519        lowLim = killBelowEnergyForA[1];
00520        
00521   if (     particle->GetDefinition() == instance->GetIon("alpha++")
00522        ||  particle->GetDefinition() == instance->GetIon("alpha+")
00523        ||  particle->GetDefinition() == instance->GetIon("helium")
00524      )
00525 
00526        lowLim = killBelowEnergyForA[4];*/
00527 
00528     G4double k = particle->GetKineticEnergy();
00529 
00530     const G4String& particleName = particle->GetDefinition()->GetParticleName();
00531 
00532     // SI - the following is useless since lowLim is already defined
00533     /*
00534   std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
00535   pos1 = lowEnergyLimit.find(particleName);
00536 
00537   if (pos1 != lowEnergyLimit.end())
00538   {
00539     lowLim = pos1->second;
00540   }
00541   */
00542 
00543     std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
00544     pos2 = highEnergyLimit.find(particleName);
00545 
00546     if (pos2 != highEnergyLimit.end())highLim = pos2->second;
00547 
00548     if (k >= lowLim && k < highLim)
00549     {
00550         G4ParticleDefinition* definition = particle->GetDefinition();
00551         G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
00552         /*
00553       G4double particleMass = definition->GetPDGMass();
00554       G4double totalEnergy = k + particleMass;
00555       G4double pSquare = k*(totalEnergy+particleMass);
00556       G4double totalMomentum = std::sqrt(pSquare);
00557       */
00558 
00559         G4int ionizationShell = RandomSelect(k,particleName);
00560 
00561         // sample deexcitation
00562         // here we assume that H_{2}O electronic levels are the same of Oxigen.
00563         // this can be considered true with a rough 10% error in energy on K-shell,
00564 
00565         G4int secNumberInit = 0;  // need to know at a certain point the enrgy of secondaries
00566         G4int secNumberFinal = 0; // So I'll make the diference and then sum the energies
00567         G4double bindingEnergy = 0;
00568         bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
00569 
00570         if(fAtomDeexcitation) {
00571             G4int Z = 8;
00572             G4AtomicShellEnumerator as = fKShell;
00573 
00574             if (ionizationShell <5 && ionizationShell >1)
00575             {
00576                 as = G4AtomicShellEnumerator(4-ionizationShell);
00577             }
00578             else if (ionizationShell <2)
00579             {
00580                 as = G4AtomicShellEnumerator(3);
00581             }
00582 
00583 
00584             //  DEBUG
00585             //  if (ionizationShell == 4) {
00586             //
00587             //    G4cout << "Z: " << Z << " as: " << as
00588             //               << " ionizationShell: " << ionizationShell << " bindingEnergy: "<< bindingEnergy/eV << G4endl;
00589             //    G4cout << "Press <Enter> key to continue..." << G4endl;
00590             //    G4cin.ignore();
00591             //  }
00592 
00593 
00594             const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
00595             secNumberInit = fvect->size();
00596             fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
00597             secNumberFinal = fvect->size();
00598         }
00599 
00600 
00601         G4double secondaryKinetic = RandomizeEjectedElectronEnergy(definition,k,ionizationShell);
00602 
00603         G4double cosTheta = 0.;
00604         G4double phi = 0.;
00605         RandomizeEjectedElectronDirection(definition, k,secondaryKinetic, cosTheta, phi, ionizationShell);
00606 
00607         G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
00608         G4double dirX = sinTheta*std::cos(phi);
00609         G4double dirY = sinTheta*std::sin(phi);
00610         G4double dirZ = cosTheta;
00611         G4ThreeVector deltaDirection(dirX,dirY,dirZ);
00612         deltaDirection.rotateUz(primaryDirection);
00613 
00614         // Ignored for ions on electrons
00615         /*
00616       G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
00617 
00618       G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
00619       G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
00620       G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
00621       G4double finalMomentum = std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz);
00622       finalPx /= finalMomentum;
00623       finalPy /= finalMomentum;
00624       finalPz /= finalMomentum;
00625 
00626       G4ThreeVector direction;
00627       direction.set(finalPx,finalPy,finalPz);
00628 
00629       fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
00630       */
00631 
00632         fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection);
00633         G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
00634         G4double deexSecEnergy = 0;
00635         for (G4int j=secNumberInit; j < secNumberFinal; j++) {
00636 
00637             deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();
00638 
00639         }
00640 
00641         fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
00642         fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy-secondaryKinetic-deexSecEnergy);
00643 
00644         G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic) ;
00645         fvect->push_back(dp);
00646 
00647         const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
00648         G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule,
00649                                                                ionizationShell,
00650                                                                theIncomingTrack);
00651     }
00652 
00653     // SI - not useful since low energy of model is 0 eV
00654 
00655     if (k < lowLim)
00656     {
00657         fParticleChangeForGamma->SetProposedKineticEnergy(0.);
00658         fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
00659         fParticleChangeForGamma->ProposeLocalEnergyDeposit(k);
00660     }
00661 
00662 }
00663 
00664 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00665 
00666 G4double G4DNARuddIonisationExtendedModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition, 
00667                                                                           G4double k,
00668                                                                           G4int shell)
00669 {
00670     //-- Fast sampling method -----
00671     G4double proposed_energy;
00672     G4double random1;
00673     G4double value_sampling;
00674     G4double max1;
00675 
00676     do
00677     {
00678         proposed_energy = ProposedSampledEnergy(particleDefinition, k, shell); // Proposed energy by inverse function sampling
00679 
00680         max1=0.;
00681 
00682         for(G4double en=0.; en<20.; en+=1.) if(RejectionFunction(particleDefinition, k, en, shell) > max1)
00683             max1=RejectionFunction(particleDefinition, k, en, shell);
00684 
00685         random1 = G4UniformRand()*max1;
00686 
00687         value_sampling = RejectionFunction(particleDefinition, k, proposed_energy, shell);
00688 
00689     } while(random1 > value_sampling);
00690 
00691     return(proposed_energy);
00692 }
00693 
00694 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00695 
00696 
00697 void G4DNARuddIonisationExtendedModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition, 
00698                                                                          G4double k,
00699                                                                          G4double secKinetic,
00700                                                                          G4double & cosTheta,
00701                                                                          G4double & phi,
00702                                                                          G4int shell )
00703 {
00704     G4double maxSecKinetic = 0.;
00705     G4double maximumEnergyTransfer = 0.;
00706 
00707     /* if (particleDefinition == G4Proton::ProtonDefinition()
00708       || particleDefinition == instance->GetIon("hydrogen"))
00709   {
00710       if(k/MeV < 100.)maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
00711       else {
00712         G4double beta2 = 1.-(1.+k*);
00713         maxSecKinetic =
00714         }
00715   }
00716   
00717   if (particleDefinition == instance->GetIon("helium")
00718       || particleDefinition == instance->GetIon("alpha+")
00719       || particleDefinition == instance->GetIon("alpha++"))
00720   {
00721       maxSecKinetic = 4.* (0.511 / 3728) * k;
00722   }
00723   
00724   if (particleDefinition == G4Carbon::Carbon())
00725   {
00726       maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k / 12.;
00727   }*/
00728 
00729     // ZF. generalized & relativistic version
00730 
00731     if( (k/MeV)/(particleDefinition->GetPDGMass()/MeV)  <= 0.1 )
00732     {
00733         maximumEnergyTransfer= 4.* (electron_mass_c2 / particleDefinition->GetPDGMass()) * k;
00734         maximumEnergyTransfer+=waterStructure.IonisationEnergy(shell);
00735     }
00736     else
00737     {
00738         G4double approx_nuc_number = particleDefinition->GetPDGMass() / proton_mass_c2;
00739         G4double en_per_nucleon = k/approx_nuc_number;
00740         G4double beta2 = 1. - 1./pow( (1.+(en_per_nucleon/electron_mass_c2)*(electron_mass_c2/proton_mass_c2)), 2.);
00741         G4double gamma = 1./sqrt(1.-beta2);
00742         maximumEnergyTransfer = 2.*electron_mass_c2*(gamma*gamma-1.)/(1.+2.*gamma*(electron_mass_c2/particleDefinition->GetPDGMass())+pow(electron_mass_c2/particleDefinition->GetPDGMass(), 2.) );
00743         maximumEnergyTransfer+=waterStructure.IonisationEnergy(shell);
00744     }
00745 
00746     maxSecKinetic = maximumEnergyTransfer-waterStructure.IonisationEnergy(shell);
00747 
00748     phi = twopi * G4UniformRand();
00749 
00750     if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
00751     else cosTheta = (2.*G4UniformRand())-1.;
00752 
00753 }
00754 
00755 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00756 
00757 G4double G4DNARuddIonisationExtendedModel::RejectionFunction(G4ParticleDefinition* particleDefinition, 
00758                                                              G4double k,
00759                                                              G4double proposed_ws,
00760                                                              G4int ionizationLevelIndex)
00761 {
00762     const G4int j=ionizationLevelIndex;
00763     G4double Bj_energy, alphaConst;
00764     G4double Ry = 13.6*eV;
00765     const G4double Gj[5] = {0.99, 1.11, 1.11, 0.52, 1.};
00766 
00767     // const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV}; //Ding Paper
00768 
00769     // Following values provided by M. Dingfelder (priv. comm)
00770     const G4double Bj[5] = {12.60*eV, 14.70*eV, 18.40*eV, 32.20*eV, 540*eV};
00771 
00772     if (j == 4)
00773     {
00774         alphaConst = 0.66;
00775         //---Note that the following (j==4) cases are provided by   M. Dingfelder (priv. comm)
00776         Bj_energy = waterStructure.IonisationEnergy(ionizationLevelIndex);
00777         //---
00778     }
00779     else
00780     {
00781         alphaConst = 0.64;
00782         Bj_energy = Bj[ionizationLevelIndex];
00783     }
00784 
00785     G4double energyTransfer = proposed_ws + Bj_energy;
00786     proposed_ws/=Bj_energy;
00787     G4DNAGenericIonsManager *instance;
00788     instance = G4DNAGenericIonsManager::Instance();
00789     G4double tau = 0.;
00790     G4double A_ion = 0.;
00791     tau = (electron_mass_c2 / particleDefinition->GetPDGMass()) * k;
00792     A_ion = particleDefinition->GetAtomicMass();
00793 
00794     G4double v2;
00795     G4double beta2;
00796 
00797     if((tau/MeV)<5.447761194e-2)
00798     {
00799         v2 = tau / Bj_energy;
00800         beta2 = 2.*tau / electron_mass_c2;
00801     }
00802     // Relativistic
00803     else
00804     {
00805         v2 = (electron_mass_c2 / 2. / Bj_energy) * (1. - (1./ pow( (1.+ (tau/electron_mass_c2)),2) ));
00806         beta2 =1. - 1./(1.+ (tau/electron_mass_c2/A_ion))/(1.+ (tau/electron_mass_c2/A_ion));
00807     }
00808 
00809     G4double v = std::sqrt(v2);
00810     G4double wc = 4.*v2 - 2.*v - (Ry/(4.*Bj_energy));
00811     G4double rejection_term = 1.+std::exp(alphaConst*(proposed_ws - wc) / v);
00812     rejection_term = (1./rejection_term)*CorrectionFactor(particleDefinition,k,ionizationLevelIndex) * Gj[j];
00813     //* (S/Bj_energy) ; Not needed anymore
00814 
00815     G4bool isHelium = false;
00816 
00817     if (    particleDefinition == G4Proton::ProtonDefinition()
00818             || particleDefinition == instance->GetIon("hydrogen")
00819             )
00820     {
00821         return(rejection_term);
00822     }
00823 
00824     else if(particleDefinition->GetAtomicMass() > 4) // anything above Helium
00825     {
00826         G4double Z = particleDefinition->GetAtomicNumber();
00827         G4double x = 100.*std::sqrt(beta2)/std::pow(Z,(2./3.));
00828         G4double Zeffion = Z*(1.-std::exp(-1.316*x+0.112*x*x-0.0650*x*x*x));
00829         rejection_term*=Zeffion*Zeffion;
00830     }
00831 
00832     else if (particleDefinition == instance->GetIon("alpha++") )
00833     {
00834         isHelium = true;
00835         slaterEffectiveCharge[0]=0.;
00836         slaterEffectiveCharge[1]=0.;
00837         slaterEffectiveCharge[2]=0.;
00838         sCoefficient[0]=0.;
00839         sCoefficient[1]=0.;
00840         sCoefficient[2]=0.;
00841     }
00842 
00843     else if (particleDefinition == instance->GetIon("alpha+") )
00844     {
00845         isHelium = true;
00846         slaterEffectiveCharge[0]=2.0;
00847         // The following values are provided by M. Dingfelder (priv. comm)
00848         slaterEffectiveCharge[1]=2.0;
00849         slaterEffectiveCharge[2]=2.0;
00850         //
00851         sCoefficient[0]=0.7;
00852         sCoefficient[1]=0.15;
00853         sCoefficient[2]=0.15;
00854     }
00855 
00856     else if (particleDefinition == instance->GetIon("helium") )
00857     {
00858         isHelium = true;
00859         slaterEffectiveCharge[0]=1.7;
00860         slaterEffectiveCharge[1]=1.15;
00861         slaterEffectiveCharge[2]=1.15;
00862         sCoefficient[0]=0.5;
00863         sCoefficient[1]=0.25;
00864         sCoefficient[2]=0.25;
00865     }
00866 
00867 //    if (    particleDefinition == instance->GetIon("helium")
00868 //            || particleDefinition == instance->GetIon("alpha+")
00869 //            || particleDefinition == instance->GetIon("alpha++")
00870 //            )
00871     if (isHelium)
00872     {
00873 
00874         G4double zEff = particleDefinition->GetPDGCharge() / eplus + particleDefinition->GetLeptonNumber();
00875 
00876         zEff -= ( sCoefficient[0] * S_1s(k, energyTransfer, slaterEffectiveCharge[0], 1.) +
00877                   sCoefficient[1] * S_2s(k, energyTransfer, slaterEffectiveCharge[1], 2.) +
00878                   sCoefficient[2] * S_2p(k, energyTransfer, slaterEffectiveCharge[2], 2.) );
00879 
00880         rejection_term*= zEff * zEff;
00881     }
00882 
00883     return (rejection_term);
00884 }
00885 
00886 
00887 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00888 
00889 
00890 G4double G4DNARuddIonisationExtendedModel::ProposedSampledEnergy(G4ParticleDefinition* particle, 
00891                                                                  G4double k,
00892                                                                  G4int ionizationLevelIndex)
00893 {
00894 
00895     const G4int j=ionizationLevelIndex;
00896 
00897     G4double A1, B1, C1, D1, E1, A2, B2, C2, D2;
00898     //G4double alphaConst ;
00899     G4double Bj_energy;
00900 
00901     // const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV}; //Ding Paper
00902     // Following values provided by M. Dingfelder (priv. comm)
00903 
00904     const G4double Bj[5] = {12.60*eV, 14.70*eV, 18.40*eV, 32.20*eV, 540*eV};
00905 
00906     if (j == 4)
00907     {
00908         //Data For Liquid Water K SHELL from Dingfelder (Protons in Water)
00909         A1 = 1.25;
00910         B1 = 0.5;
00911         C1 = 1.00;
00912         D1 = 1.00;
00913         E1 = 3.00;
00914         A2 = 1.10;
00915         B2 = 1.30;
00916         C2 = 1.00;
00917         D2 = 0.00;
00918         //alphaConst = 0.66;
00919         //---Note that the following (j==4) cases are provided by   M. Dingfelder (priv. comm)
00920         Bj_energy = waterStructure.IonisationEnergy(ionizationLevelIndex);
00921         //---
00922     }
00923     else
00924     {
00925         //Data For Liquid Water from Dingfelder (Protons in Water)
00926         A1 = 1.02;
00927         B1 = 82.0;
00928         C1 = 0.45;
00929         D1 = -0.80;
00930         E1 = 0.38;
00931         A2 = 1.07;
00932         //B2 = 14.6; From Ding Paper
00933         // Value provided by M. Dingfelder (priv. comm)
00934         B2 = 11.6;
00935         //
00936         C2 = 0.60;
00937         D2 = 0.04;
00938         //alphaConst = 0.64;
00939 
00940         Bj_energy = Bj[ionizationLevelIndex];
00941     }
00942 
00943     G4double tau = 0.;
00944     G4double A_ion = 0.;
00945     tau = (electron_mass_c2 / particle->GetPDGMass()) * k;
00946 
00947     A_ion = particle->GetAtomicMass();
00948 
00949     G4double v2;
00950     G4double beta2;
00951     if((tau/MeV)<5.447761194e-2)
00952     {
00953         v2 = tau / Bj_energy;
00954         beta2 = 2.*tau / electron_mass_c2;
00955     }
00956     // Relativistic
00957     else
00958     {
00959         v2 = (electron_mass_c2 / 2. / Bj_energy) * (1. - (1./ pow( (1.+ (tau/electron_mass_c2)),2) ));
00960         beta2 =1. - 1./(1.+ (tau/electron_mass_c2/A_ion))/(1.+ (tau/electron_mass_c2/A_ion));
00961     }
00962 
00963     G4double v = std::sqrt(v2);
00964     //G4double wc = 4.*v2 - 2.*v - (Ry/(4.*Bj_energy));
00965     G4double L1 = (C1* std::pow(v,(D1))) / (1.+ E1*std::pow(v, (D1+4.)));
00966     G4double L2 = C2*std::pow(v,(D2));
00967     G4double H1 = (A1*std::log(1.+v2)) / (v2+(B1/v2));
00968     G4double H2 = (A2/v2) + (B2/(v2*v2));
00969     G4double F1 = L1+H1;
00970     G4double F2 = (L2*H2)/(L2+H2);
00971 
00972     // ZF. generalized & relativistic version
00973     G4double maximumEnergy;
00974 
00975     //---- maximum kinetic energy , non relativistic ------
00976     if( (k/MeV)/(particle->GetPDGMass()/MeV)  <= 0.1 )
00977     {
00978         maximumEnergy = 4.* (electron_mass_c2 / particle->GetPDGMass()) * k;
00979     }
00980     //---- relativistic -----------------------------------
00981     else
00982     {
00983         G4double gamma = 1./sqrt(1.-beta2);
00984         maximumEnergy = 2.*electron_mass_c2*(gamma*gamma-1.)/
00985                 (1.+2.*gamma*(electron_mass_c2/particle->GetPDGMass())+pow(electron_mass_c2/particle->GetPDGMass(), 2.) );
00986     }
00987 
00988     //either it is transfered energy or secondary electron energy ...
00989     //maximumEnergy-=Bj_energy;
00990 
00991     //-----------------------------------------------------
00992     G4double wmax = maximumEnergy/Bj_energy;
00993     G4double c = wmax*(F2*wmax+F1*(2.+wmax))/(2.*(1.+wmax)*(1.+wmax));
00994     c=1./c; 
00995     G4double randVal = G4UniformRand();
00996     G4double proposed_ws = F1*F1*c*c + 2.*F2*c*randVal - 2.*F1*c*randVal;
00997     proposed_ws = -F1*c+2.*randVal+std::sqrt(proposed_ws);
00998     //  proposed_ws = -F1*c+2.*randVal-std::sqrt(proposed_ws);
00999     proposed_ws/= ( F1*c + F2*c - 2.*randVal );
01000     proposed_ws*=Bj_energy;
01001 
01002     return(proposed_ws);
01003 
01004 }
01005 
01006 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
01007 
01008 G4double G4DNARuddIonisationExtendedModel::S_1s(G4double t, 
01009                                                 G4double energyTransferred,
01010                                                 G4double slaterEffectiveChg,
01011                                                 G4double shellNumber)
01012 {
01013     // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2)
01014     // Dingfelder, in Chattanooga 2005 proceedings, formula (7)
01015 
01016     G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
01017     G4double value = 1. - std::exp(-2 * r) * ( ( 2. * r + 2. ) * r + 1. );
01018 
01019     return value;
01020 }
01021 
01022 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
01023 
01024 G4double G4DNARuddIonisationExtendedModel::S_2s(G4double t,
01025                                                 G4double energyTransferred,
01026                                                 G4double slaterEffectiveChg,
01027                                                 G4double shellNumber)
01028 {
01029     // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4)
01030     // Dingfelder, in Chattanooga 2005 proceedings, formula (8)
01031 
01032     G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
01033     G4double value =  1. - std::exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
01034 
01035     return value;
01036 
01037 }
01038 
01039 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
01040 
01041 G4double G4DNARuddIonisationExtendedModel::S_2p(G4double t, 
01042                                                 G4double energyTransferred,
01043                                                 G4double slaterEffectiveChg,
01044                                                 G4double shellNumber)
01045 {
01046     // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^3 + 2/3 r^4)
01047     // Dingfelder, in Chattanooga 2005 proceedings, formula (9)
01048 
01049     G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
01050     G4double value =  1. - std::exp(-2 * r) * (((( 2./3. * r + 4./3.) * r + 2.) * r + 2.) * r  + 1.);
01051 
01052     return value;
01053 }
01054 
01055 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
01056 
01057 G4double G4DNARuddIonisationExtendedModel::R(G4double t,
01058                                              G4double energyTransferred,
01059                                              G4double slaterEffectiveChg,
01060                                              G4double shellNumber)
01061 {
01062     // tElectron = m_electron / m_alpha * t
01063     // Dingfelder, in Chattanooga 2005 proceedings, p 4
01064 
01065     G4double tElectron = 0.511/3728. * t;
01066     // The following values are provided by M. Dingfelder (priv. comm)
01067     G4double H = 2.*13.60569172 * eV;
01068     G4double value = std::sqrt ( 2. * tElectron / H ) / ( energyTransferred / H ) *  (slaterEffectiveChg/shellNumber);
01069 
01070     return value;
01071 }
01072 
01073 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
01074 
01075 G4double G4DNARuddIonisationExtendedModel::CorrectionFactor(G4ParticleDefinition* particleDefinition, G4double k, G4int shell) 
01076 {
01077     // ZF Shortened
01078     G4DNAGenericIonsManager *instance;
01079     instance = G4DNAGenericIonsManager::Instance();
01080 
01081     if (particleDefinition == instance->GetIon("hydrogen") && shell < 4)
01082     {
01083         G4double value = (std::log10(k/eV)-4.2)/0.5;
01084         // The following values are provided by M. Dingfelder (priv. comm)
01085         return((0.6/(1+std::exp(value))) + 0.9);
01086     }
01087     else
01088     {
01089         return(1.);
01090     }
01091 }
01092 
01093 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
01094 
01095 G4int G4DNARuddIonisationExtendedModel::RandomSelect(G4double k, const G4String& particle )
01096 {   
01097 
01098     G4int level = 0;
01099 
01100     // Retrieve data table corresponding to the current particle type
01101 
01102     std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
01103     pos = tableData.find(particle);
01104 
01105     if (pos != tableData.end())
01106     {
01107         G4DNACrossSectionDataSet* table = pos->second;
01108 
01109         if (table != 0)
01110         {
01111             G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
01112 
01113             const size_t n(table->NumberOfComponents());
01114             size_t i(n);
01115             G4double value = 0.;
01116 
01117             while (i>0)
01118             {
01119                 i--;
01120                 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
01121 
01122                 value += valuesBuffer[i];
01123             }
01124 
01125             value *= G4UniformRand();
01126 
01127             i = n;
01128 
01129             while (i > 0)
01130             {
01131                 i--;
01132 
01133                 if (valuesBuffer[i] > value)
01134                 {
01135                     delete[] valuesBuffer;
01136                     return i;
01137                 }
01138                 value -= valuesBuffer[i];
01139             }
01140 
01141             if (valuesBuffer) delete[] valuesBuffer;
01142 
01143         }
01144     }
01145     else
01146     {
01147         G4Exception("G4DNARuddIonisationExtendedModel::RandomSelect","em0002",
01148                     FatalException,"Model not applicable to particle type.");
01149     }
01150 
01151     return level;
01152 }
01153 
01154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
01155 
01156 G4double G4DNARuddIonisationExtendedModel::PartialCrossSection(const G4Track& track )
01157 {
01158     G4double sigma = 0.;
01159 
01160     const G4DynamicParticle* particle = track.GetDynamicParticle();
01161     G4double k = particle->GetKineticEnergy();
01162 
01163     G4double lowLim = 0;
01164     G4double highLim = 0;
01165 
01166     const G4String& particleName = particle->GetDefinition()->GetParticleName();
01167 
01168     std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
01169     pos1 = lowEnergyLimit.find(particleName);
01170 
01171     if (pos1 != lowEnergyLimit.end())
01172     {
01173         lowLim = pos1->second;
01174     }
01175 
01176     std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
01177     pos2 = highEnergyLimit.find(particleName);
01178 
01179     if (pos2 != highEnergyLimit.end())
01180     {
01181         highLim = pos2->second;
01182     }
01183 
01184     if (k >= lowLim && k <= highLim)
01185     {
01186         std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
01187         pos = tableData.find(particleName);
01188 
01189         if (pos != tableData.end())
01190         {
01191             G4DNACrossSectionDataSet* table = pos->second;
01192             if (table != 0)
01193             {
01194                 sigma = table->FindValue(k);
01195             }
01196         }
01197         else
01198         {
01199             G4Exception("G4DNARuddIonisationExtendedModel::PartialCrossSection","em0002",
01200                         FatalException,"Model not applicable to particle type.");
01201         }
01202     }
01203 
01204     return sigma;
01205 }
01206 
01207 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
01208 
01209 G4double G4DNARuddIonisationExtendedModel::Sum(G4double /* energy */, const G4String& /* particle */)
01210 {
01211     return 0;
01212 }
01213 

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