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/G4CrossSectionIonisationBornPartial.cc

    r819 r961  
    2424// ********************************************************************
    2525//
    26 //
    27 // $Id: G4CrossSectionIonisationBornPartial.cc,v 1.3 2007/11/09 20:11:04 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: G4CrossSectionIonisationBornPartial.cc,v 1.5 2009/01/20 07:40:53 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    5228
    5329#include "G4CrossSectionIonisationBornPartial.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 
    61 #include "Randomize.hh"
    62 
    63 #include <utility>
     30
     31//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    6432
    6533G4CrossSectionIonisationBornPartial::G4CrossSectionIonisationBornPartial()
    6634{
    67 
    68   name = "IonisationBorn";
    69  
    70   // Default energy limits (defined for protection against anomalous behaviour only)
    71   name = "IonisationBornPartial";
    72   lowEnergyLimitDefault = 25 * eV;
     35  lowEnergyLimitDefault = 12.61 * eV;
    7336  highEnergyLimitDefault = 30 * keV;
    7437
    … …  
    8245  G4String proton;
    8346 
    84   // Factor to scale microscopic/macroscopic cross section data in water
    85   // ---- MGP ---- Hardcoded (taken from prototype code); to be replaced with proper calculation
    8647  G4double scaleFactor = (1.e-22 / 3.343) * m*m;
    8748
    88 
    89   // Data members for electrons
    90 
    9149  if (electronDef != 0)
    92     {
    93       electron = electronDef->GetParticleName();
    94       tableFile[electron] = fileElectron;
    95 
    96       // Energy limits
    97       lowEnergyLimit[electron] = 25. * eV;
    98       highEnergyLimit[electron] = 30. * keV;
    99 
    100       // Create data set with electron cross section data and load values stored in file
    101       G4DNACrossSectionDataSet* tableE = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
    102       tableE->LoadData(fileElectron);
     50  {
     51    electron = electronDef->GetParticleName();
     52    tableFile[electron] = fileElectron;
     53
     54    lowEnergyLimit[electron] = 12.61 * eV;
     55    highEnergyLimit[electron] = 30. * keV;
     56
     57    G4DNACrossSectionDataSet* tableE = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
     58    tableE->LoadData(fileElectron);
    10359     
    104       // Insert key-table pair in map
    105       tableData[electron] = tableE;
    106     }
     60    tableData[electron] = tableE;
     61  }
    10762  else
    108     {
    109       G4Exception("G4CrossSectionIonisationBornPartial Constructor: electron is not defined");
    110     }
    111 
    112   // Data members for protons
     63  {
     64    G4Exception("G4CrossSectionIonisationBornPartial Constructor: electron is not defined");
     65  }
    11366
    11467  if (protonDef != 0)
    115     {
    116       proton = protonDef->GetParticleName();
    117       tableFile[proton] = fileProton;
    118 
    119       // Energy limits
    120       lowEnergyLimit[proton] = 500. * keV;
    121       highEnergyLimit[proton] = 10. * MeV;
    122 
    123       // Create data set with proton cross section data and load values stored in file
    124       G4DNACrossSectionDataSet* tableP = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
    125       tableP->LoadData(fileProton);
     68  {
     69    proton = protonDef->GetParticleName();
     70    tableFile[proton] = fileProton;
     71
     72    lowEnergyLimit[proton] = 500. * keV;
     73    highEnergyLimit[proton] = 10. * MeV;
     74
     75    G4DNACrossSectionDataSet* tableP = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
     76    tableP->LoadData(fileProton);
    12677     
    127       // Insert key-table pair in map
    128       tableData[proton] = tableP;
    129     }
     78    tableData[proton] = tableP;
     79  }
    13080  else
    131     {
    132       G4Exception("G4CrossSectionIonisationBornPartial Constructor: proton is not defined");
    133     }
    134 }
    135 
     81  {
     82    G4Exception("G4CrossSectionIonisationBornPartial Constructor: proton is not defined");
     83  }
     84}
     85
     86//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    13687
    13788G4CrossSectionIonisationBornPartial::~G4CrossSectionIonisationBornPartial()
    13889{
    139    // Destroy the content of the map
    14090  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
    14191  for (pos = tableData.begin(); pos != tableData.end(); ++pos)
    142     {
    143       G4DNACrossSectionDataSet* table = pos->second;
    144       delete table;
    145     }
    146 }
    147 
     92  {
     93    G4DNACrossSectionDataSet* table = pos->second;
     94    delete table;
     95  }
     96}
     97
     98//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    14899
    149100G4int G4CrossSectionIonisationBornPartial::RandomSelect(G4double k, const G4String& particle )
    150101{   
    151  
    152102  G4int level = 0;
    153103
    154   // Retrieve data table corresponding to the current particle type 
    155 
    156     std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
    157     pos = tableData.find(particle);
    158 
    159     if (pos != tableData.end())
     104  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
     105  pos = tableData.find(particle);
     106
     107  if (pos != tableData.end())
     108  {
     109    G4DNACrossSectionDataSet* table = pos->second;
     110
     111    if (table != 0)
     112    {
     113      G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
     114      const size_t n(table->NumberOfComponents());
     115      size_t i(n);
     116      G4double value = 0.;
     117           
     118      while (i>0)
     119      {
     120        i--;
     121        valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
     122        value += valuesBuffer[i];
     123      }
     124           
     125      value *= G4UniformRand();
     126   
     127      i = n;
     128           
     129      while (i > 0)
    160130      {
    161         G4DNACrossSectionDataSet* table = pos->second;
    162 
    163         if (table != 0)
    164           {
    165             // C-style arrays are used in G4DNACrossSectionDataSet: this design feature was
    166             // introduced without authorization and should be replaced by the use of STL containers
    167            
    168             G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
    169            
    170             const size_t n(table->NumberOfComponents());
    171             size_t i(n);
    172             G4double value = 0.;
    173            
    174             while (i>0)
    175               {
    176                 i--;
    177                 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
    178                 value += valuesBuffer[i];
    179               }
    180            
    181             value *= G4UniformRand();
    182            
    183             i = n;
    184            
    185             while (i > 0)
    186               {
    187                 i--;
    188                
    189                 if (valuesBuffer[i] > value)
    190                   {
    191                     delete[] valuesBuffer;
    192                     return i;
    193                   }
    194                 value -= valuesBuffer[i];
    195               }
    196            
    197             // It should never end up here
    198 
    199             // ---- MGP ---- Is the following line really necessary? 
    200             if (valuesBuffer) delete[] valuesBuffer;
    201            
    202           }
     131        i--;
     132
     133        if (valuesBuffer[i] > value)
     134        {
     135          delete[] valuesBuffer;
     136          return i;
     137        }
     138        value -= valuesBuffer[i];
    203139      }
    204     else
    205       {
    206         G4Exception("G4CrossSectionIonisationBornPartial: attempting to calculate cross section for wrong particle");
    207       }
     140           
     141      if (valuesBuffer) delete[] valuesBuffer;
     142   
     143    }
     144  }
     145  else
     146  {
     147    G4Exception("G4CrossSectionIonisationBornPartial: attempting to calculate cross section for wrong particle");
     148  }
    208149     
    209150  return level;
    210151}
    211152
     153//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    212154
    213155G4double G4CrossSectionIonisationBornPartial::CrossSection(const G4Track& track )
    … …  
    218160  G4double k = particle->GetKineticEnergy();
    219161 
    220   // Cross section = 0 outside the energy validity limits set in the constructor
    221   // ---- MGP ---- Better handling of these limits to be set in a following design iteration
    222 
    223162  G4double lowLim = lowEnergyLimitDefault;
    224163  G4double highLim = highEnergyLimitDefault;
    … …  
    226165  const G4String& particleName = particle->GetDefinition()->GetParticleName();
    227166
    228   // Retrieve energy limits for the current particle type
    229 
    230     std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
    231     pos1 = lowEnergyLimit.find(particleName);
    232 
    233     // Lower limit
    234     if (pos1 != lowEnergyLimit.end())
    235       {
    236         lowLim = pos1->second;
    237       }
    238 
    239     // Upper limit
    240     std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
    241     pos2 = highEnergyLimit.find(particleName);
    242 
    243     if (pos2 != highEnergyLimit.end())
    244       {
    245         highLim = pos2->second;
    246       }
    247 
    248     // Verify that the current track is within the energy limits of validity of the cross section model
    249     if (k > lowLim && k < highLim)
    250       {
    251         std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
    252         pos = tableData.find(particleName);
     167  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
     168  pos1 = lowEnergyLimit.find(particleName);
     169
     170  if (pos1 != lowEnergyLimit.end())
     171  {
     172    lowLim = pos1->second;
     173  }
     174
     175  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
     176  pos2 = highEnergyLimit.find(particleName);
     177
     178  if (pos2 != highEnergyLimit.end())
     179  {
     180     highLim = pos2->second;
     181  }
     182
     183  if (k > lowLim && k < highLim)
     184  {
     185     std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
     186     pos = tableData.find(particleName);
    253187       
    254         if (pos != tableData.end())
    255           {
    256             G4DNACrossSectionDataSet* table = pos->second;
    257             if (table != 0)
    258               {
    259                 // ---- MGP ---- Temporary
    260                 // table->PrintData();
    261 
    262                 // Cross section
    263                 sigma = table->FindValue(k);
    264               }
    265           }
    266         else
    267           {
    268             // The track corresponds to a not pertinent particle
    269             G4Exception("G4CrossSectionIonisationBornPartial: attempting to calculate cross section for wrong particle");
    270           }
    271       }
    272 
    273     return sigma;
    274 }
    275 
     188     if (pos != tableData.end())
     189     {
     190       G4DNACrossSectionDataSet* table = pos->second;
     191       if (table != 0)
     192       {
     193         sigma = table->FindValue(k);
     194       }
     195     }
     196     else
     197     {
     198       G4Exception("G4CrossSectionIonisationBornPartial: attempting to calculate cross section for wrong particle");
     199     }
     200  }
     201
     202  return sigma;
     203}
     204
     205//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    276206
    277207G4double G4CrossSectionIonisationBornPartial::Sum(G4double /* energy */, const G4String& /* particle */)
    278208{
    279 
    280209  return 0;
    281210}
Note: See TracChangeset for help on using the changeset viewer.