Ignore:
Timestamp:
Nov 5, 2010, 3:45:55 PM (14 years ago)
Author:
garnier
Message:

update ti head

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4DNARuddIonisationModel.cc,v 1.17 2010/04/07 20:08:31 sincerti Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4DNARuddIonisationModel.cc,v 1.21 2010/11/04 14:52:17 sincerti Exp $
     27// GEANT4 tag $Name: emlowen-V09-03-54 $
    2828//
    2929
    … …  
    159159
    160160    lowEnergyLimit[alphaPlusPlus] = lowEnergyLimitForZ2;
    161     highEnergyLimit[alphaPlusPlus] = 10. * MeV;
     161    highEnergyLimit[alphaPlusPlus] = 400. * MeV;
    162162
    163163    // Cross section
    … …  
    181181
    182182    lowEnergyLimit[alphaPlus] = lowEnergyLimitForZ2;
    183     highEnergyLimit[alphaPlus] = 10. * MeV;
     183    highEnergyLimit[alphaPlus] = 400. * MeV;
    184184
    185185    // Cross section
    … …  
    202202
    203203    lowEnergyLimit[helium] = lowEnergyLimitForZ2;
    204     highEnergyLimit[helium] = 10. * MeV;
     204    highEnergyLimit[helium] = 400. * MeV;
    205205
    206206    // Cross section
    … …  
    363363              sigma = table->FindValue(k);
    364364
    365               // BEGIN ELECTRON CORRECTION
    366               // add ONE or TWO electron-water excitation for alpha+ and helium
    367    
    368               if ( particleDefinition == instance->GetIon("alpha+")
    369                    ||
    370                    particleDefinition == instance->GetIon("helium")
    371                    )
    372               {
    373      
    374                   G4DNACrossSectionDataSet* electronDataset = new G4DNACrossSectionDataSet
    375                     (new G4LogLogInterpolation, eV, (1./3.343e22)*m*m);
    376        
    377                   electronDataset->LoadData("dna/sigma_ionisation_e_born");
    378 
    379                   G4double kElectron = k * 0.511/3728;
    380 
    381                   if ( particleDefinition == instance->GetIon("alpha+") )
    382                   {
    383                       G4double tmp1 = table->FindValue(k) + electronDataset->FindValue(kElectron);
    384                       delete electronDataset;
    385                       if (verboseLevel > 3)
    386                       {
    387                         G4cout << "---> Kinetic energy(eV)=" << k/eV << G4endl;
    388                         G4cout << " - Cross section per water molecule (cm^2)=" << tmp1/cm/cm << G4endl;
    389                         G4cout << " - Cross section per water molecule (cm^-1)=" <<
    390                         tmp1*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
    391                       }
    392                       return tmp1*material->GetAtomicNumDensityVector()[1];
    393                   }
    394 
    395                   if ( particleDefinition == instance->GetIon("helium") )
    396                   {
    397                       G4double tmp2 = table->FindValue(k) +  2. * electronDataset->FindValue(kElectron);
    398                       delete electronDataset;
    399                       if (verboseLevel > 3)
    400                       {
    401                         G4cout << "---> Kinetic energy(eV)=" << k/eV << G4endl;
    402                         G4cout << " - Cross section per water molecule (cm^2)=" << tmp2/cm/cm << G4endl;
    403                         G4cout << " - Cross section per water molecule (cm^-1)=" << tmp2*
    404                         material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
    405                       }
    406                       return tmp2*material->GetAtomicNumDensityVector()[1];
    407                   }
    408               }     
    409 
    410               // END ELECTRON CORRECTION
    411365         }
    412366      }
    … …  
    489443      G4ParticleDefinition* definition = particle->GetDefinition();
    490444      G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
     445      /*
    491446      G4double particleMass = definition->GetPDGMass();
    492447      G4double totalEnergy = k + particleMass;
    493448      G4double pSquare = k*(totalEnergy+particleMass);
    494449      G4double totalMomentum = std::sqrt(pSquare);
     450      */
    495451
    496452      G4int ionizationShell = RandomSelect(k,particleName);
    … …  
    511467      deltaDirection.rotateUz(primaryDirection);
    512468
     469      // Ignored for ions on electrons
     470      /*
    513471      G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
    514472
    … …  
    525483
    526484      fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
     485      */
     486      fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection);
     487
    527488      fParticleChangeForGamma->SetProposedKineticEnergy(k-bindingEnergy-secondaryKinetic);
    528489      fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy);
    … …  
    575536  }
    576537 
     538 
    577539  G4double secElecKinetic = 0.;
    578540 
    … …  
    616578 
    617579  phi = twopi * G4UniformRand();
    618   cosTheta = std::sqrt(secKinetic / maxSecKinetic);
     580 
     581  //cosTheta = std::sqrt(secKinetic / maxSecKinetic);
     582
     583  // Restriction below 100 eV from Emfietzoglou (2000)
     584
     585  if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
     586  else cosTheta = (2.*G4UniformRand())-1.;
     587 
    619588}
    620589
    … …  
    656625  G4double alphaConst ;
    657626
    658   const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV};
     627//  const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV};
     628// The following values are provided by M. dingfelder (priv. comm)
     629  const G4double Bj[5] = {12.60*eV, 14.70*eV, 18.40*eV, 32.20*eV, 540*eV};
    659630
    660631  if (j == 4)
    … …  
    681652      E1 = 0.38;
    682653      A2 = 1.07;
    683       B2 = 14.6;
     654      // Value provided by M. Dingfelder (priv. comm)
     655      B2 = 11.6;
     656      //
    684657      C2 = 0.60;
    685658      D2 = 0.04;
    … …  
    697670
    698671  G4double w = wBig / Bj[ionizationLevelIndex];
     672  // Note that the following (j==4) cases are provided by   M. Dingfelder (priv. comm)
     673  if (j==4) w = wBig / waterStructure.IonisationEnergy(ionizationLevelIndex);
     674
    699675  G4double Ry = 13.6*eV;
    700676
    … …  
    713689      tau = (0.511/3728.) * k ;
    714690  }
    715  
     691
    716692  G4double S = 4.*pi * Bohr_radius*Bohr_radius * n * std::pow((Ry/Bj[ionizationLevelIndex]),2);
     693  if (j==4) S = 4.*pi * Bohr_radius*Bohr_radius * n * std::pow((Ry/waterStructure.IonisationEnergy(ionizationLevelIndex)),2);
     694
    717695  G4double v2 = tau / Bj[ionizationLevelIndex];
     696  if (j==4) v2 = tau / waterStructure.IonisationEnergy(ionizationLevelIndex);
     697
    718698  G4double v = std::sqrt(v2);
    719699  G4double wc = 4.*v2 - 2.*v - (Ry/(4.*Bj[ionizationLevelIndex]));
     700  if (j==4) wc = 4.*v2 - 2.*v - (Ry/(4.*waterStructure.IonisationEnergy(ionizationLevelIndex)));
    720701
    721702  G4double L1 = (C1* std::pow(v,(D1))) / (1.+ E1*std::pow(v, (D1+4.)));
    … …  
    731712    * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
    732713
     714  if (j==4) sigma = CorrectionFactor(particleDefinition, k)
     715    * Gj[j] * (S/waterStructure.IonisationEnergy(ionizationLevelIndex))
     716    * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
     717
    733718  if ( (particleDefinition == instance->GetIon("hydrogen")) && (ionizationLevelIndex==4))
    734719
    735     sigma = Gj[j] * (S/Bj[ionizationLevelIndex])
     720//    sigma = Gj[j] * (S/Bj[ionizationLevelIndex])
     721    sigma = Gj[j] * (S/waterStructure.IonisationEnergy(ionizationLevelIndex))
    736722    * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
    737723
    … …  
    756742  {
    757743      slaterEffectiveCharge[0]=2.0;
    758       slaterEffectiveCharge[1]=1.15;
    759       slaterEffectiveCharge[2]=1.15;
     744      // The following values are provided by M. Dingfelder (priv. comm)   
     745      slaterEffectiveCharge[1]=2.0;
     746      slaterEffectiveCharge[2]=2.0;
     747      //
    760748      sCoefficient[0]=0.7;
    761749      sCoefficient[1]=0.15;
    … …  
    780768      sigma = Gj[j] * (S/Bj[ionizationLevelIndex]) * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
    781769   
     770      if (j==4) sigma = Gj[j] * (S/waterStructure.IonisationEnergy(ionizationLevelIndex))
     771                              * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
     772
    782773      G4double zEff = particleDefinition->GetPDGCharge() / eplus + particleDefinition->GetLeptonNumber();
    783774 
    … …  
    852843
    853844  G4double tElectron = 0.511/3728. * t;
    854   G4double value = 2. * tElectron * slaterEffectiveChg / (energyTransferred * shellNumber);
     845  // The following values are provided by M. Dingfelder (priv. comm)   
     846  G4double H = 2.*13.60569172 * eV;
     847  G4double value = std::sqrt ( 2. * tElectron / H ) / ( energyTransferred / H ) *  (slaterEffectiveChg/shellNumber);
    855848 
    856849  return value;
    … …  
    872865    {
    873866        G4double value = (std::log10(k/eV)-4.2)/0.5;
    874         return((0.8/(1+std::exp(value))) + 0.9);
     867        // The following values are provided by M. Dingfelder (priv. comm)   
     868        return((0.6/(1+std::exp(value))) + 0.9);
    875869    }
    876870    else
    … …  
    891885  G4DNAGenericIonsManager *instance;
    892886  instance = G4DNAGenericIonsManager::Instance();
    893   G4double kElectron(0);
    894 
    895   G4DNACrossSectionDataSet * electronDataset = new G4DNACrossSectionDataSet (new G4LogLogInterpolation, eV, (1./3.343e22)*m*m);
    896887 
    897   if ( particle == instance->GetIon("alpha+")->GetParticleName()
    898        ||
    899        particle == instance->GetIon("helium")->GetParticleName()
    900        )
    901   {     
    902       electronDataset->LoadData("dna/sigma_ionisation_e_born");
    903 
    904       kElectron = k * 0.511/3728;
    905        
    906   }     
    907  
    908   // END PART 1/2 OF ELECTRON CORRECTION
    909  
    910888  G4int level = 0;
    911889
    … …  
    931909              i--;
    932910              valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
    933 
    934               // BEGIN PART 2/2 OF ELECTRON CORRECTION
    935               // Use only electron partial cross sections
    936              
    937               if (particle == instance->GetIon("alpha+")->GetParticleName())
    938                 {valuesBuffer[i]=table->GetComponent(i)->FindValue(k) + electronDataset->GetComponent(i)->FindValue(kElectron); }
    939 
    940               if (particle == instance->GetIon("helium")->GetParticleName())
    941                 {valuesBuffer[i]=table->GetComponent(i)->FindValue(k) + 2*electronDataset->GetComponent(i)->FindValue(kElectron); }
    942 
    943               // BEGIN PART 2/2 OF ELECTRON CORRECTION
    944 
    945911              value += valuesBuffer[i];
    946912          }
    … …  
    958924              {
    959925                  delete[] valuesBuffer;
    960                  
    961                   if (electronDataset) delete electronDataset;
    962                  
    963926                  return i;
    964927              }
    … …  
    975938  }
    976939
    977   delete electronDataset;
    978      
    979940  return level;
    980941}
Note: See TracChangeset for help on using the changeset viewer.