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

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