Ignore:
Timestamp:
Apr 6, 2009, 12:21:12 PM (15 years ago)
Author:
garnier
Message:

update processes

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/electromagnetic/lowenergy/src/G4CrossSectionIonisationRudd.cc

    r819 r961  
    2424// ********************************************************************
    2525//
    26 //
    27 // $Id: G4CrossSectionIonisationRudd.cc,v 1.2 2007/11/09 16:20:16 pia Exp $
    28 // GEANT4 tag $Name:  $
    29 //
    30 // Contact Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
    31 //
    32 // Reference: TNS Geant4-DNA paper
    33 // Reference for implementation model: NIM. 155, pp. 145-156, 1978
    34 
    35 // History:
    36 // -----------
    37 // Date         Name              Modification
    38 // 28 Apr 2007  M.G. Pia          Created in compliance with design described in TNS paper
    39 //
    40 // -------------------------------------------------------------------
    41 
    42 // Class description:
    43 // Geant4-DNA Cross total cross section for electron elastic scattering in water
    44 // Reference: TNS Geant4-DNA paper
    45 // S. Chauvie et al., Geant4 physics processes for microdosimetry simulation:
    46 // design foundation and implementation of the first set of models,
    47 // IEEE Trans. Nucl. Sci., vol. 54, no. 6, Dec. 2007.
    48 // Further documentation available from http://www.ge.infn.it/geant4/dna
    49 
    50 // -------------------------------------------------------------------
    51 
     26// $Id: G4CrossSectionIonisationRudd.cc,v 1.4 2008/07/14 20:47:34 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    5228
    5329#include "G4CrossSectionIonisationRudd.hh"
    54 #include "G4ParticleDefinition.hh"
    55 #include "G4Electron.hh"
    56 #include "G4Proton.hh"
    57 #include "G4Track.hh"
    58 #include "G4LogLogInterpolation.hh"
    59 #include "G4SystemOfUnits.hh"
    60 #include "G4DNAGenericIonsManager.hh"
     30
     31//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    6132
    6233G4CrossSectionIonisationRudd::G4CrossSectionIonisationRudd()
    6334{
    64 
    65   name = "IonisationRudd";
    66  
    67   // Default energy limits (defined for protection against anomalous behaviour only)
    68   // ZERO LOW ENERGY LIMIT FOR ALLOWED KILLING
    6935  lowEnergyLimitDefault = 0 * eV;
    7036  highEnergyLimitDefault = 100 * MeV;
    … …  
    9056  G4String helium;
    9157
    92   // Factor to scale microscopic/macroscopic cross section data in water
    93   // ---- MGP ---- Hardcoded (taken from prototype code); to be replaced with proper calculation
    9458  G4double scaleFactor = 1 * m*m;
    9559
    96   // Data members for protons
    97 
    9860  if (protonDef != 0)
    99     {
    100       proton = protonDef->GetParticleName();
    101       tableFile[proton] = fileProton;
    102 
    103       // Energy limits
    104       lowEnergyLimit[proton] = 0. * eV;
    105       highEnergyLimit[proton] = 500. * keV;
    106 
    107       // Create data set with proton cross section data and load values stored in file
    108       G4DNACrossSectionDataSet* tableProton = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
     61  {
     62    proton = protonDef->GetParticleName();
     63    tableFile[proton] = fileProton;
     64
     65    lowEnergyLimit[proton] = 0. * eV;
     66    highEnergyLimit[proton] = 500. * keV;
     67
     68    G4DNACrossSectionDataSet* tableProton = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
    10969                                                                           eV,
    11070                                                                           scaleFactor );
    111       tableProton->LoadData(fileProton);
    112      
    113       // Insert key-table pair in map
    114       tableData[proton] = tableProton;
    115     }
    116   else
    117     {
    118       G4Exception("G4CrossSectionIonisationRudd Constructor: proton is not defined");
    119     }
    120 
    121   // Data members for hydrogen
     71    tableProton->LoadData(fileProton);
     72    tableData[proton] = tableProton;
     73  }
     74  else
     75  {
     76    G4Exception("G4CrossSectionIonisationRudd Constructor: proton is not defined");
     77  }
    12278
    12379  if (hydrogenDef != 0)
    124     {
    125       hydrogen = hydrogenDef->GetParticleName();
    126       tableFile[hydrogen] = fileHydrogen;
    127 
    128       // Energy limits
    129       lowEnergyLimit[hydrogen] = 0. * eV;
    130       highEnergyLimit[hydrogen] = 100. * MeV;
    131 
    132       // Create data set with hydrogen cross section data and load values stored in file
    133       G4DNACrossSectionDataSet* tableHydrogen = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
     80  {
     81    hydrogen = hydrogenDef->GetParticleName();
     82    tableFile[hydrogen] = fileHydrogen;
     83
     84    lowEnergyLimit[hydrogen] = 0. * eV;
     85    highEnergyLimit[hydrogen] = 100. * MeV;
     86
     87    G4DNACrossSectionDataSet* tableHydrogen = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
    13488                                                                             eV,
    13589                                                                             scaleFactor );
    136       tableHydrogen->LoadData(fileHydrogen);
     90    tableHydrogen->LoadData(fileHydrogen);
    13791     
    138       // Insert key-table pair in map
    139       tableData[hydrogen] = tableHydrogen;
    140     }
    141   else
    142     {
    143       G4Exception("G4CrossSectionIonisationRudd Constructor: hydrogen is not defined");
    144     }
    145 
    146   // Data members for alphaPlusPlus
     92    tableData[hydrogen] = tableHydrogen;
     93  }
     94  else
     95  {
     96    G4Exception("G4CrossSectionIonisationRudd Constructor: hydrogen is not defined");
     97  }
    14798
    14899  if (alphaPlusPlusDef != 0)
    149     {
    150       alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
    151       tableFile[alphaPlusPlus] = fileAlphaPlusPlus;
    152 
    153       // Energy limits
    154       lowEnergyLimit[alphaPlusPlus] = 0. * keV;
    155       highEnergyLimit[alphaPlusPlus] = 10. * MeV;
    156 
    157       // Create data set with hydrogen cross section data and load values stored in file
    158       G4DNACrossSectionDataSet* tableAlphaPlusPlus = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
     100  {
     101    alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
     102    tableFile[alphaPlusPlus] = fileAlphaPlusPlus;
     103
     104    lowEnergyLimit[alphaPlusPlus] = 0. * eV;
     105    highEnergyLimit[alphaPlusPlus] = 10. * MeV;
     106
     107    G4DNACrossSectionDataSet* tableAlphaPlusPlus = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
    159108                                                                                  eV,
    160109                                                                                  scaleFactor );
    161       tableAlphaPlusPlus->LoadData(fileAlphaPlusPlus);
     110    tableAlphaPlusPlus->LoadData(fileAlphaPlusPlus);
    162111     
    163       // Insert key-table pair in map
    164       tableData[alphaPlusPlus] = tableAlphaPlusPlus;
    165     }
    166   else
    167     {
    168       G4Exception("G4CrossSectionIonisationRudd Constructor: alphaPlusPlus is not defined");
    169     }
    170 
    171   // Data members for alphaPlus
     112    tableData[alphaPlusPlus] = tableAlphaPlusPlus;
     113  }
     114  else
     115  {
     116    G4Exception("G4CrossSectionIonisationRudd Constructor: alphaPlusPlus is not defined");
     117  }
    172118
    173119  if (alphaPlusDef != 0)
    174     {
    175       alphaPlus = alphaPlusDef->GetParticleName();
    176       tableFile[alphaPlus] = fileAlphaPlus;
    177 
    178       // Energy limits
    179       lowEnergyLimit[alphaPlus] = 0. * keV;
    180       highEnergyLimit[alphaPlus] = 10. * MeV;
    181 
    182       // Create data set with hydrogen cross section data and load values stored in file
    183       G4DNACrossSectionDataSet* tableAlphaPlus = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
     120  {
     121    alphaPlus = alphaPlusDef->GetParticleName();
     122    tableFile[alphaPlus] = fileAlphaPlus;
     123
     124    lowEnergyLimit[alphaPlus] = 0. * eV;
     125    highEnergyLimit[alphaPlus] = 10. * MeV;
     126
     127    G4DNACrossSectionDataSet* tableAlphaPlus = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
    184128                                                                              eV,
    185129                                                                              scaleFactor );
    186       tableAlphaPlus->LoadData(fileAlphaPlus);
    187      
    188       // Insert key-table pair in map
    189       tableData[alphaPlus] = tableAlphaPlus;
    190     }
    191   else
    192     {
    193       G4Exception("G4CrossSectionIonisationRudd Constructor: alphaPlus is not defined");
    194     }
    195 
    196   // Data members for helium
     130    tableAlphaPlus->LoadData(fileAlphaPlus);
     131    tableData[alphaPlus] = tableAlphaPlus;
     132  }
     133  else
     134  {
     135    G4Exception("G4CrossSectionIonisationRudd Constructor: alphaPlus is not defined");
     136  }
    197137
    198138  if (heliumDef != 0)
    199     {
    200       helium = heliumDef->GetParticleName();
    201       tableFile[helium] = fileHelium;
    202 
    203       // Energy limits
    204       lowEnergyLimit[helium] = 0. * keV;
    205       highEnergyLimit[helium] = 10. * MeV;
    206 
    207       // Create data set with hydrogen cross section data and load values stored in file
    208       G4DNACrossSectionDataSet* tableHelium = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
     139  {
     140    helium = heliumDef->GetParticleName();
     141    tableFile[helium] = fileHelium;
     142
     143    lowEnergyLimit[helium] = 0. * eV;
     144    highEnergyLimit[helium] = 10. * MeV;
     145
     146    G4DNACrossSectionDataSet* tableHelium = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
    209147                                                                           eV,
    210148                                                                           scaleFactor );
    211       tableHelium->LoadData(fileHelium);
    212      
    213       // Insert key-table pair in map
    214       tableData[helium] = tableHelium;
     149    tableHelium->LoadData(fileHelium);
     150    tableData[helium] = tableHelium;
    215151    }
    216152  else
    217     {
    218       G4Exception("G4CrossSectionIonisationRudd Constructor: helium is not defined");
    219     }
     153  {
     154    G4Exception("G4CrossSectionIonisationRudd Constructor: helium is not defined");
     155  }
    220156}
    221157
     158//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    222159
    223160G4CrossSectionIonisationRudd::~G4CrossSectionIonisationRudd()
    224161{
    225   // Destroy the content of the map
    226162  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
    227163  for (pos = tableData.begin(); pos != tableData.end(); ++pos)
    228     {
    229       G4DNACrossSectionDataSet* table = pos->second;
    230       delete table;
    231     }
     164  {
     165    G4DNACrossSectionDataSet* table = pos->second;
     166    delete table;
     167  }
    232168}
    233169
    234 
     170//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    235171
    236172G4double G4CrossSectionIonisationRudd::CrossSection(const G4Track& track )
    … …  
    252188      &&
    253189      track.GetDefinition() != instance->GetIon("helium")
    254       )
     190     )
    255191           
    256192    G4Exception("G4CrossSectionIonisationRudd: attempting to calculate cross section for wrong particle");
    257193
    258   // Cross section = 0 outside the energy validity limits set in the constructor
    259194  G4double sigma = 0.;
    260195
    261   // ---- MGP ---- Better handling of these limits to be set in a following design iteration
    262196  G4double lowLim = lowEnergyLimitDefault;
    263197  G4double highLim = highEnergyLimitDefault;
    … …  
    265199  const G4String& particleName = particle->GetDefinition()->GetParticleName();
    266200
    267   // Retrieve energy limits for the current particle type
    268 
    269201  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
    270202  pos1 = lowEnergyLimit.find(particleName);
    271203
    272   // Lower limit
    273204  if (pos1 != lowEnergyLimit.end())
    274     {
    275       lowLim = pos1->second;
    276     }
    277 
    278   // Upper limit
     205  {
     206    lowLim = pos1->second;
     207  }
     208
    279209  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
    280210  pos2 = highEnergyLimit.find(particleName);
    281211
    282212  if (pos2 != highEnergyLimit.end())
     213  {
     214    highLim = pos2->second;
     215  }
     216
     217  if (k >= lowLim && k <= highLim)
     218  {
     219    std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
     220    pos = tableData.find(particleName);
     221       
     222    if (pos != tableData.end())
    283223    {
    284       highLim = pos2->second;
    285     }
    286 
    287   // Verify that the current track is within the energy limits of validity of the cross section model
    288   if (k >= lowLim && k <= highLim)
    289     {
    290       std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
    291       pos = tableData.find(particleName);
    292        
    293       if (pos != tableData.end())
    294         {
    295           G4DNACrossSectionDataSet* table = pos->second;
    296           if (table != 0)
    297             {
    298               // Cross section
     224       G4DNACrossSectionDataSet* table = pos->second;
     225       if (table != 0)
     226       {
    299227              sigma = table->FindValue(k);
    300 
    301228
    302229              // BEGIN ELECTRON CORRECTION
    … …  
    307234                   particle->GetDefinition() == instance->GetIon("helium")
    308235                   )
    309                 {
     236              {
    310237     
    311238                  G4DNACrossSectionDataSet* electronDataset = new G4DNACrossSectionDataSet
    … …  
    317244
    318245                  if ( particle->GetDefinition() == instance->GetIon("alpha+") )
    319                     {
     246                  {
    320247                      G4double tmp1 = table->FindValue(k) + electronDataset->FindValue(kElectron);
    321248                      delete electronDataset;
    322249                      return tmp1;
    323                     }
     250                  }
    324251
    325252                  if ( particle->GetDefinition() == instance->GetIon("helium") )
    326                     {
     253                  {
    327254                      G4double tmp2 = table->FindValue(k) +  2. * electronDataset->FindValue(kElectron);
    328255                      delete electronDataset;
    329256                      return tmp2;
    330                     }
    331                 }     
     257                  }
     258              }     
    332259
    333260              // END ELECTRON CORRECTION
    334 
    335             }
    336         }
    337       else
    338         {
    339           // The track does not corresponds to a particle pertinent to this model
    340           G4Exception("G4CrossSectionIonisationRudd: attempting to calculate cross section for wrong particle");
    341         }
     261       }
    342262    }
     263    else
     264    {
     265        G4Exception("G4CrossSectionIonisationRudd: attempting to calculate cross section for wrong particle");
     266    }
     267  }
    343268
    344269  return sigma;
Note: See TracChangeset for help on using the changeset viewer.