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

    r819 r961  
    2424// ********************************************************************
    2525//
    26 //
    27 // $Id: G4FinalStateIonisationBorn.cc,v 1.9 2007/11/26 17:27:09 pia Exp $
    28 // GEANT4 tag $Name:  $
    29 //
    30 // Contact Author: Sebastien Incerti (incerti@cenbg.in2p3.fr)
    31 //                 Maria Grazia Pia  (Maria.Grazia.Pia@cern.ch)
    32 //
    33 // Reference: TNS Geant4-DNA paper
    34 // Reference for implementation model: NIM. 155, pp. 145-156, 1978
    35 
    36 // History:
    37 // -----------
    38 // Date         Name              Modification
    39 // 28 Apr 2007  M.G. Pia          Created in compliance with design described in TNS paper
    40 //    Nov 2007  S. Incerti        Implementation
    41 // 26 Nov 2007  MGP               Cleaned up std::
    42 //
    43 // -------------------------------------------------------------------
    44 
    45 // Class description:
    46 // Reference: TNS Geant4-DNA paper
    47 // S. Chauvie et al., Geant4 physics processes for microdosimetry simulation:
    48 // design foundation and implementation of the first set of models,
    49 // IEEE Trans. Nucl. Sci., vol. 54, no. 6, Dec. 2007.
    50 // Further documentation available from http://www.ge.infn.it/geant4/dna
    51 
    52 // -------------------------------------------------------------------
    53 
     26// $Id: G4FinalStateIonisationBorn.cc,v 1.16 2008/12/06 13:47:12 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    5428
    5529#include "G4FinalStateIonisationBorn.hh"
    56 #include "G4Track.hh"
    57 #include "G4Step.hh"
    58 #include "G4DynamicParticle.hh"
    59 #include "Randomize.hh"
    60 
    61 #include "G4ParticleTypes.hh"
    62 #include "G4ParticleDefinition.hh"
    63 #include "G4Electron.hh"
    64 #include "G4Proton.hh"
    65 #include "G4SystemOfUnits.hh"
    66 #include "G4ParticleMomentum.hh"
    67 
     30
     31//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    6832
    6933G4FinalStateIonisationBorn::G4FinalStateIonisationBorn()
    7034{
    71 
    72   name = "IonisationBorn";
    73 
    74   // NEW
    75   // Factor to scale microscopic/macroscopic cross section data in water
    76  
    7735  G4double scaleFactor = (1.e-22 / 3.343) * m*m;
    7836
    79   // Energy limits
    8037  G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
    8138  G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
    … …  
    8441  G4String proton;
    8542
    86   // Default energy limits (defined for protection against anomalous behaviour only)
    87   lowEnergyLimitDefault = 25 * eV;
     43  lowEnergyLimitDefault = 12.61 * eV; // SI: i/o 25 eV
    8844  highEnergyLimitDefault = 10 * MeV;
    8945
    … …  
    9349    G4Exception("G4DNACrossSectionDataSet::FullFileName - G4LEDATA environment variable not set");
    9450
    95   // Data members for electrons
    96 
    9751  if (electronDef != 0)
    98     {
    99       electron = electronDef->GetParticleName();
    100       lowEnergyLimit[electron] = 25. * eV;
    101       highEnergyLimit[electron] = 30. * keV;
    102 
    103       std::ostringstream eFullFileName;
    104       eFullFileName << path << "/dna/sigmadiff_ionisation_e_born.dat";
    105       std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
    106       // eDiffCrossSection(eFullFileName.str().c_str());
    107       if (!eDiffCrossSection)
    108         {
    109           // G4cout << "ERROR OPENING DATA FILE IN ELECTRON BORN IONIZATION !!! " << G4endl;
    110           G4Exception("G4FinalStateIonisationBorn::ERROR OPENING electron DATA FILE");
    111           while(1); // ---- MGP ---- What is this?
    112         }
     52  {
     53    electron = electronDef->GetParticleName();
     54    lowEnergyLimit[electron] = 12.61 * eV; // SI: i/o 25 eV
     55    highEnergyLimit[electron] = 30. * keV;
     56
     57    std::ostringstream eFullFileName;
     58    eFullFileName << path << "/dna/sigmadiff_ionisation_e_born.dat";
     59    std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
     60    if (!eDiffCrossSection)
     61    {
     62      G4Exception("G4FinalStateIonisationBorn::ERROR OPENING electron DATA FILE");
     63    }
    11364     
    114       eTdummyVec.push_back(0.);
    115       while(!eDiffCrossSection.eof())
    116         {
    117           double tDummy;
    118           double eDummy;
    119           eDiffCrossSection>>tDummy>>eDummy;
    120           if (tDummy != eTdummyVec.back()) eTdummyVec.push_back(tDummy);
    121           for (int j=0; j<5; j++)
    122             {
    123               eDiffCrossSection>>eDiffCrossSectionData[j][tDummy][eDummy];
    124               eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
    125               eVecm[tDummy].push_back(eDummy);
    126             }
    127         }
    128 
    129     }
     65    eTdummyVec.push_back(0.);
     66    while(!eDiffCrossSection.eof())
     67    {
     68      double tDummy;
     69      double eDummy;
     70      eDiffCrossSection>>tDummy>>eDummy;
     71      if (tDummy != eTdummyVec.back()) eTdummyVec.push_back(tDummy);
     72      for (int j=0; j<5; j++)
     73      {
     74        eDiffCrossSection>>eDiffCrossSectionData[j][tDummy][eDummy];
     75
     76        // SI - only if eof is not reached !
     77        if (!eDiffCrossSection.eof()) eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
     78
     79        eVecm[tDummy].push_back(eDummy);
     80
     81      }
     82    }
     83
     84  }
    13085  else
    131     {
    132       G4Exception("G4FinalStateIonisationBorn Constructor: electron is not defined");
    133     }
    134 
    135   // Data members for protons
     86  {
     87    G4Exception("G4FinalStateIonisationBorn Constructor: electron is not defined");
     88  }
    13689
    13790  if (protonDef != 0)
    138     {
    139       proton = protonDef->GetParticleName();
    140       lowEnergyLimit[proton] = 500. * keV;
    141       highEnergyLimit[proton] = 10. * MeV;
    142 
    143       std::ostringstream pFullFileName;
    144       pFullFileName << path << "/dna/sigmadiff_ionisation_p_born.dat";
    145       std::ifstream pDiffCrossSection(pFullFileName.str().c_str());
    146       //     pDiffCrossSection(pFullFileName.str().c_str());
    147       if (!pDiffCrossSection)
    148         {
    149           // G4cout<<"ERROR OPENING DATA FILE IN PROTON BORN IONIZATION !!! "<<G4endl;
    150           G4Exception("G4FinalStateIonisationBorn::ERROR OPENING proton DATA FILE");
    151           while(1); // ---- MGP ---- What is this?
    152         }
     91  {
     92    proton = protonDef->GetParticleName();
     93    lowEnergyLimit[proton] = 500. * keV;
     94    highEnergyLimit[proton] = 10. * MeV;
     95
     96    std::ostringstream pFullFileName;
     97    pFullFileName << path << "/dna/sigmadiff_ionisation_p_born.dat";
     98    std::ifstream pDiffCrossSection(pFullFileName.str().c_str());
     99    if (!pDiffCrossSection)
     100    {
     101      G4Exception("G4FinalStateIonisationBorn::ERROR OPENING proton DATA FILE");
     102    }
    153103     
    154       pTdummyVec.push_back(0.);
    155       while(!pDiffCrossSection.eof())
    156         {
    157           double tDummy;
    158           double eDummy;
    159           pDiffCrossSection>>tDummy>>eDummy;
    160           if (tDummy != pTdummyVec.back()) pTdummyVec.push_back(tDummy);
    161           for (int j=0; j<5; j++)
    162             {
    163               pDiffCrossSection>>pDiffCrossSectionData[j][tDummy][eDummy];
    164               pDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
    165               //G4cout << "j=" << j << " Tdum=" << tDummy << " Edum=" << eDummy << " pDiff=" << pDiffCrossSectionData[j][tDummy][eDummy] << G4endl;
    166               pVecm[tDummy].push_back(eDummy);
    167             }
    168         }
    169     }
     104    pTdummyVec.push_back(0.);
     105    while(!pDiffCrossSection.eof())
     106    {
     107      double tDummy;
     108      double eDummy;
     109      pDiffCrossSection>>tDummy>>eDummy;
     110      if (tDummy != pTdummyVec.back()) pTdummyVec.push_back(tDummy);
     111      for (int j=0; j<5; j++)
     112      {
     113        pDiffCrossSection>>pDiffCrossSectionData[j][tDummy][eDummy];
     114
     115        // SI - only if eof is not reached !
     116        if (!pDiffCrossSection.eof()) pDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
     117
     118        pVecm[tDummy].push_back(eDummy);
     119      }
     120    }
     121  }
    170122  else
    171     {
    172       G4Exception("G4FinalStateIonisationBorn Constructor: proton is not defined");
    173     }
    174 }
    175 
     123  {
     124    G4Exception("G4FinalStateIonisationBorn Constructor: proton is not defined");
     125  }
     126}
     127
     128//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    176129
    177130G4FinalStateIonisationBorn::~G4FinalStateIonisationBorn()
    … …  
    181134}
    182135
     136//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    183137
    184138const G4FinalStateProduct& G4FinalStateIonisationBorn::GenerateFinalState(const G4Track& track, const G4Step& /* step */)
    185139{
    186   // Clear previous secondaries, energy deposit and particle kill status
    187140  product.Clear();
    188141
    … …  
    196149  const G4String& particleName = particle->GetDefinition()->GetParticleName();
    197150
    198   // Retrieve energy limits for the current particle type
    199 
    200151  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
    201152  pos1 = lowEnergyLimit.find(particleName);
    202153
    203   // Lower limit
    204154  if (pos1 != lowEnergyLimit.end())
    205     {
    206       lowLim = pos1->second;
    207     }
    208 
    209   // Upper limit
     155  {
     156    lowLim = pos1->second;
     157  }
     158
    210159  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
    211160  pos2 = highEnergyLimit.find(particleName);
    212161
    213162  if (pos2 != highEnergyLimit.end())
    214     {
    215       highLim = pos2->second;
    216     }
    217 
    218   // Verify that the current track is within the energy limits of validity of the cross section model
     163  {
     164    highLim = pos2->second;
     165  }
    219166
    220167  if (k >= lowLim && k <= highLim)
    221     {
    222       // Kinetic energy of primary particle
    223 
    224       G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
    225       G4double particleMass = particle->GetDefinition()->GetPDGMass();
    226       G4double totalEnergy = k + particleMass;
    227       G4double pSquare = k * (totalEnergy + particleMass);
    228       G4double totalMomentum = std::sqrt(pSquare);
    229 
    230       const G4String& particleName = particle->GetDefinition()->GetParticleName();
    231  
    232       G4int ionizationShell = cross.RandomSelect(k,particleName);
    233  
    234       G4double secondaryKinetic = RandomizeEjectedElectronEnergy(particle->GetDefinition(),k,ionizationShell);
    235  
    236       G4double bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
    237 
    238       G4double cosTheta = 0.;
    239       G4double phi = 0.;
    240       RandomizeEjectedElectronDirection(track.GetDefinition(), k,secondaryKinetic, cosTheta, phi);
    241 
    242       G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
    243       G4double dirX = sinTheta*std::cos(phi);
    244       G4double dirY = sinTheta*std::sin(phi);
    245       G4double dirZ = cosTheta;
    246       G4ThreeVector deltaDirection(dirX,dirY,dirZ);
    247       deltaDirection.rotateUz(primaryDirection);
    248 
    249       G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
    250 
    251       //Primary Particle Direction
    252       G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
    253       G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
    254       G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
    255       G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
    256       finalPx /= finalMomentum;
    257       finalPy /= finalMomentum;
    258       finalPz /= finalMomentum;
    259 
    260       product.ModifyPrimaryParticle(finalPx,finalPy,finalPz,k-bindingEnergy-secondaryKinetic);
    261       product.AddEnergyDeposit(bindingEnergy);
    262 
    263       G4DynamicParticle* aElectron = new G4DynamicParticle(G4Electron::Electron(),deltaDirection,secondaryKinetic);
    264       product.AddSecondary(aElectron);
    265     }
     168  {
     169    G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
     170    G4double particleMass = particle->GetDefinition()->GetPDGMass();
     171    G4double totalEnergy = k + particleMass;
     172    G4double pSquare = k * (totalEnergy + particleMass);
     173    G4double totalMomentum = std::sqrt(pSquare);
     174
     175    const G4String& particleName = particle->GetDefinition()->GetParticleName();
     176 
     177    G4int ionizationShell = cross.RandomSelect(k,particleName);
     178 
     179    G4double secondaryKinetic = RandomizeEjectedElectronEnergy(particle->GetDefinition(),k,ionizationShell);
     180 
     181    G4double bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
     182
     183    G4double cosTheta = 0.;
     184    G4double phi = 0.;
     185    RandomizeEjectedElectronDirection(track.GetDefinition(), k,secondaryKinetic, cosTheta, phi);
     186
     187    G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
     188    G4double dirX = sinTheta*std::cos(phi);
     189    G4double dirY = sinTheta*std::sin(phi);
     190    G4double dirZ = cosTheta;
     191    G4ThreeVector deltaDirection(dirX,dirY,dirZ);
     192    deltaDirection.rotateUz(primaryDirection);
     193
     194    G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
     195
     196    G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
     197    G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
     198    G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
     199    G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
     200    finalPx /= finalMomentum;
     201    finalPy /= finalMomentum;
     202    finalPz /= finalMomentum;
     203
     204    product.ModifyPrimaryParticle(finalPx,finalPy,finalPz,k-bindingEnergy-secondaryKinetic);
     205    product.AddEnergyDeposit(bindingEnergy);
     206
     207    G4DynamicParticle* aElectron = new G4DynamicParticle(G4Electron::Electron(),deltaDirection,secondaryKinetic);
     208    product.AddSecondary(aElectron);
     209  }
    266210
    267211  return product;
    268212}
    269213
     214//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    270215
    271216G4double G4FinalStateIonisationBorn::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
    272                                                                     G4double k,
    273                                                                     G4int shell)
    274 {
    275 
     217G4double k, G4int shell)
     218{
    276219  if (particleDefinition == G4Electron::ElectronDefinition())
    277     {
    278 
    279       G4double maximumEnergyTransfer=0.;
    280       if ((k+waterStructure.IonisationEnergy(shell))/2. > k) maximumEnergyTransfer=k;
    281       else maximumEnergyTransfer = (k+waterStructure.IonisationEnergy(shell))/2.;
     220  {
     221    G4double maximumEnergyTransfer=0.;
     222    if ((k+waterStructure.IonisationEnergy(shell))/2. > k) maximumEnergyTransfer=k;
     223    else maximumEnergyTransfer = (k+waterStructure.IonisationEnergy(shell))/2.;
    282224   
    283       G4double crossSectionMaximum = 0.;
    284       for(G4double value=waterStructure.IonisationEnergy(shell); value<=maximumEnergyTransfer; value+=0.1*eV)
    285         {
    286           G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
    287           if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
    288         }
     225    G4double crossSectionMaximum = 0.;
     226    for(G4double value=waterStructure.IonisationEnergy(shell); value<=maximumEnergyTransfer; value+=0.1*eV)
     227    {
     228      G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
     229      if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
     230    }
    289231 
    290       G4double secondaryElectronKineticEnergy=0.;
    291       do
    292         {
    293           secondaryElectronKineticEnergy = G4UniformRand() * (maximumEnergyTransfer-waterStructure.IonisationEnergy(shell));
    294         } while(G4UniformRand()*crossSectionMaximum >
     232    G4double secondaryElectronKineticEnergy=0.;
     233    do
     234    {
     235      secondaryElectronKineticEnergy = G4UniformRand() * (maximumEnergyTransfer-waterStructure.IonisationEnergy(shell));
     236    } while(G4UniformRand()*crossSectionMaximum >
     237      DifferentialCrossSection(particleDefinition, k/eV,(secondaryElectronKineticEnergy+waterStructure.IonisationEnergy(shell))/eV,shell));
     238
     239    return secondaryElectronKineticEnergy;
     240 
     241  }
     242 
     243  if (particleDefinition == G4Proton::ProtonDefinition())
     244  {
     245    G4double maximumKineticEnergyTransfer = 4.* (electron_mass_c2 / proton_mass_c2) * k - (waterStructure.IonisationEnergy(shell));
     246
     247    G4double crossSectionMaximum = 0.;
     248    for (G4double value = waterStructure.IonisationEnergy(shell);
     249         value<=4.*waterStructure.IonisationEnergy(shell) ;
     250         value+=0.1*eV)
     251    {
     252      G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
     253      if (differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
     254    }
     255
     256    G4double secondaryElectronKineticEnergy = 0.;
     257    do
     258    {
     259      secondaryElectronKineticEnergy = G4UniformRand() * maximumKineticEnergyTransfer;
     260    } while(G4UniformRand()*crossSectionMaximum >=
    295261              DifferentialCrossSection(particleDefinition, k/eV,(secondaryElectronKineticEnergy+waterStructure.IonisationEnergy(shell))/eV,shell));
    296262
    297       return secondaryElectronKineticEnergy;
    298  
    299     }
    300  
    301   if (particleDefinition == G4Proton::ProtonDefinition())
    302     {
    303       G4double maximumKineticEnergyTransfer = 4.* (electron_mass_c2 / proton_mass_c2) * k - (waterStructure.IonisationEnergy(shell));
    304 
    305       G4double crossSectionMaximum = 0.;
    306       for (G4double value = waterStructure.IonisationEnergy(shell);
    307            value<=4.*waterStructure.IonisationEnergy(shell) ;
    308            value+=0.1*eV)
    309         {
    310           G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
    311           if (differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
    312         }
    313 
    314       G4double secondaryElectronKineticEnergy = 0.;
    315       do
    316         {
    317           secondaryElectronKineticEnergy = G4UniformRand() * maximumKineticEnergyTransfer;
    318         } while(G4UniformRand()*crossSectionMaximum >=
    319               DifferentialCrossSection(particleDefinition, k/eV,(secondaryElectronKineticEnergy+waterStructure.IonisationEnergy(shell))/eV,shell));
    320 
    321       return secondaryElectronKineticEnergy;
    322     }
     263    return secondaryElectronKineticEnergy;
     264  }
    323265
    324266  return 0;
    325267}
    326268
     269//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    327270
    328271void G4FinalStateIonisationBorn::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
    … …  
    333276{
    334277  if (particleDefinition == G4Electron::ElectronDefinition())
    335     {
    336 
    337       phi = twopi * G4UniformRand();
    338       if (secKinetic < 50.*eV) cosTheta = (2.*G4UniformRand())-1.;
    339       else if (secKinetic <= 200.*eV)   
    340         {
    341           if (G4UniformRand() <= 0.1) cosTheta = (2.*G4UniformRand())-1.;
    342           else cosTheta = G4UniformRand()*(std::sqrt(2.)/2);
    343         }
    344       else     
    345         {
    346           G4double sin2O = (1.-secKinetic/k) / (1.+secKinetic/(2.*electron_mass_c2));
    347           cosTheta = std::sqrt(1.-sin2O);
    348         }
    349     }
     278  {
     279    phi = twopi * G4UniformRand();
     280    if (secKinetic < 50.*eV) cosTheta = (2.*G4UniformRand())-1.;
     281    else if (secKinetic <= 200.*eV)     
     282    {
     283      if (G4UniformRand() <= 0.1) cosTheta = (2.*G4UniformRand())-1.;
     284      else cosTheta = G4UniformRand()*(std::sqrt(2.)/2);
     285    }
     286    else       
     287    {
     288      G4double sin2O = (1.-secKinetic/k) / (1.+secKinetic/(2.*electron_mass_c2));
     289      cosTheta = std::sqrt(1.-sin2O);
     290    }
     291  }
    350292 
    351293  if (particleDefinition == G4Proton::ProtonDefinition())
    352     {
    353       G4double maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
    354       phi = twopi * G4UniformRand();
    355       cosTheta = std::sqrt(secKinetic / maxSecKinetic);
    356     }                   
    357 }
    358 
     294  {
     295    G4double maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
     296    phi = twopi * G4UniformRand();
     297    cosTheta = std::sqrt(secKinetic / maxSecKinetic);
     298  }                     
     299}
     300
     301//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    359302
    360303double G4FinalStateIonisationBorn::DifferentialCrossSection(G4ParticleDefinition * particleDefinition,
    … …  
    366309
    367310  if (energyTransfer >= waterStructure.IonisationEnergy(ionizationLevelIndex))
    368     {
    369       G4double valueT1 = 0;
    370       G4double valueT2 = 0;
    371       G4double valueE21 = 0;
    372       G4double valueE22 = 0;
    373       G4double valueE12 = 0;
    374       G4double valueE11 = 0;
    375 
    376       G4double xs11  =  0;   
    377       G4double xs12 = 0;
    378       G4double xs21 = 0;
    379       G4double xs22 = 0;
    380 
     311  {
     312    G4double valueT1 = 0;
     313    G4double valueT2 = 0;
     314    G4double valueE21 = 0;
     315    G4double valueE22 = 0;
     316    G4double valueE12 = 0;
     317    G4double valueE11 = 0;
     318
     319    G4double xs11 = 0;   
     320    G4double xs12 = 0;
     321    G4double xs21 = 0;
     322    G4double xs22 = 0;
     323
     324    if (particleDefinition == G4Electron::ElectronDefinition())
     325    {
     326      // k should be in eV and energy transfer eV also
     327
     328      std::vector<double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
     329
     330      std::vector<double>::iterator t1 = t2-1;
     331
     332      // SI : the following condition avoids situations where energyTransfer >last vector element
     333      if (energyTransfer <= eVecm[(*t1)].back())
     334      {
     335        std::vector<double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),eVecm[(*t1)].end(), energyTransfer);
     336        std::vector<double>::iterator e11 = e12-1;
     337
     338        std::vector<double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),eVecm[(*t2)].end(), energyTransfer);
     339        std::vector<double>::iterator e21 = e22-1;
     340
     341        valueT1  =*t1;
     342        valueT2  =*t2;
     343        valueE21 =*e21;
     344        valueE22 =*e22;
     345        valueE12 =*e12;
     346        valueE11 =*e11;
     347
     348        xs11 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
     349        xs12 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
     350        xs21 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
     351        xs22 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
     352      }
     353
     354    }
     355 
     356   if (particleDefinition == G4Proton::ProtonDefinition())
     357   {
     358      // k should be in eV and energy transfer eV also
     359      std::vector<double>::iterator t2 = std::upper_bound(pTdummyVec.begin(),pTdummyVec.end(), k);
     360      std::vector<double>::iterator t1 = t2-1;
     361     
     362        std::vector<double>::iterator e12 = std::upper_bound(pVecm[(*t1)].begin(),pVecm[(*t1)].end(), energyTransfer);
     363        std::vector<double>::iterator e11 = e12-1;
     364
     365        std::vector<double>::iterator e22 = std::upper_bound(pVecm[(*t2)].begin(),pVecm[(*t2)].end(), energyTransfer);
     366        std::vector<double>::iterator e21 = e22-1;
    381367 
    382       if (particleDefinition == G4Electron::ElectronDefinition())
    383         {
    384           // k should be in eV and energy transfer eV also
    385           std::vector<double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
    386           std::vector<double>::iterator t1 = t2-1;
    387           std::vector<double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),eVecm[(*t1)].end(), energyTransfer);
    388           std::vector<double>::iterator e11 = e12-1;
    389 
    390           std::vector<double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),eVecm[(*t2)].end(), energyTransfer);
    391           std::vector<double>::iterator e21 = e22-1;
    392 
    393           valueT1  =*t1;
    394           valueT2  =*t2;
    395           valueE21 =*e21;
    396           valueE22 =*e22;
    397           valueE12 =*e12;
    398           valueE11 =*e11;
    399 
    400           xs11 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
    401           xs12 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
    402           xs21 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
    403           xs22 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
    404 
    405         }
    406  
    407       if (particleDefinition == G4Proton::ProtonDefinition())
    408         {
    409           // k should be in eV and energy transfer eV also
    410           std::vector<double>::iterator t2 = std::upper_bound(pTdummyVec.begin(),pTdummyVec.end(), k);
    411           std::vector<double>::iterator t1 = t2-1;
    412           std::vector<double>::iterator e12 = std::upper_bound(pVecm[(*t1)].begin(),pVecm[(*t1)].end(), energyTransfer);
    413           std::vector<double>::iterator e11 = e12-1;
    414 
    415           std::vector<double>::iterator e22 = std::upper_bound(pVecm[(*t2)].begin(),pVecm[(*t2)].end(), energyTransfer);
    416           std::vector<double>::iterator e21 = e22-1;
    417  
    418           valueT1  =*t1;
    419           valueT2  =*t2;
    420           valueE21 =*e21;
    421           valueE22 =*e22;
    422           valueE12 =*e12;
    423           valueE11 =*e11;
    424 
    425           xs11 = pDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
    426           xs12 = pDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
    427           xs21 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
    428           xs22 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
    429         }
    430  
    431       G4double xsProduct = xs11 * xs12 * xs21 * xs22;
    432       // if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) return (0.);
    433       if (xsProduct != 0.)
    434         {
    435           sigma = QuadInterpolator(valueE11, valueE12,
     368        valueT1  =*t1;
     369        valueT2  =*t2;
     370        valueE21 =*e21;
     371        valueE22 =*e22;
     372        valueE12 =*e12;
     373        valueE11 =*e11;
     374
     375        xs11 = pDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
     376        xs12 = pDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
     377        xs21 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
     378        xs22 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
     379
     380   }
     381 
     382   G4double xsProduct = xs11 * xs12 * xs21 * xs22;
     383   if (xsProduct != 0.)
     384   {
     385     sigma = QuadInterpolator(     valueE11, valueE12,
    436386                                   valueE21, valueE22,
    437387                                   xs11, xs12,
    … …  
    439389                                   valueT1, valueT2,
    440390                                   k, energyTransfer);
    441         }
    442     }
     391   }
     392 
     393 }
     394 
    443395  return sigma;
    444396}
    445397
     398//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    446399
    447400G4double G4FinalStateIonisationBorn::LogLogInterpolate(G4double e1,
    … …  
    458411}
    459412
     413//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    460414
    461415G4double G4FinalStateIonisationBorn::QuadInterpolator(G4double e11, G4double e12,
Note: See TracChangeset for help on using the changeset viewer.