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/standard/src/G4eCoulombScatteringModel.cc

    r819 r961  
    2424// ********************************************************************
    2525//
    26 // $Id: G4eCoulombScatteringModel.cc,v 1.40 2008/01/07 08:32:01 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: G4eCoulombScatteringModel.cc,v 1.59 2008/10/22 18:39:29 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    4444// 19.08.06 V.Ivanchenko add inline function ScreeningParameter
    4545// 09.10.07 V.Ivanchenko reorganized methods, add cut dependence in scattering off e-
     46// 09.06.08 V.Ivanchenko add SelectIsotope and sampling of the recoil ion
    4647//
    4748// Class Description:
    … …  
    6162#include "G4Positron.hh"
    6263#include "G4Proton.hh"
     64#include "G4ParticleTable.hh"
    6365
    6466//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    6668using namespace std;
    6769
    68 G4eCoulombScatteringModel::G4eCoulombScatteringModel(
    69   G4double thetaMin, G4double thetaMax, G4bool build,
    70   G4double tlim, const G4String& nam)
     70G4eCoulombScatteringModel::G4eCoulombScatteringModel(const G4String& nam)
    7171  : G4VEmModel(nam),
    72     cosThetaMin(cos(thetaMin)),
    73     cosThetaMax(cos(thetaMax)),
    74     q2Limit(tlim),
    75     theCrossSectionTable(0),
    76     lowKEnergy(keV),
    77     highKEnergy(TeV),
     72    cosThetaMin(1.0),
     73    cosThetaMax(-1.0),
     74    q2Limit(TeV*TeV),
    7875    alpha2(fine_structure_const*fine_structure_const),
    7976    faclim(100.0),
    80     nbins(12),
    81     nmax(100),
    82     buildTable(build),
    8377    isInitialised(false)
    8478{
    8579  fNistManager = G4NistManager::Instance();
     80  theParticleTable = G4ParticleTable::GetParticleTable();
    8681  theElectron = G4Electron::Electron();
    8782  thePositron = G4Positron::Positron();
    8883  theProton   = G4Proton::Proton();
     84  currentMaterial = 0;
     85  currentElement  = 0;
    8986  a0 = alpha2*electron_mass_c2*electron_mass_c2/(0.885*0.885);
    9087  G4double p0 = electron_mass_c2*classic_electr_radius;
    9188  coeff  = twopi*p0*p0;
    9289  constn = 6.937e-6/(MeV*MeV);
    93   tkin = targetZ = targetA = mom2 = DBL_MIN;
     90  tkin = targetZ = mom2 = DBL_MIN;
    9491  elecXSection = nucXSection = 0.0;
     92  recoilThreshold = DBL_MAX;
    9593  ecut = DBL_MAX;
    9694  particle = 0;
    97   for(size_t j=0; j<100; j++) {index[j] = -1;}
     95  currentCouple = 0;
     96  for(size_t j=0; j<100; j++) {
     97    FF[j] = 0.0;
     98  }
    9899}
    99100
    … …  
    101102
    102103G4eCoulombScatteringModel::~G4eCoulombScatteringModel()
    103 {
    104   if(theCrossSectionTable) {
    105     theCrossSectionTable->clearAndDestroy();
    106     delete theCrossSectionTable;
    107   }
    108 }
     104{}
    109105
    110106//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    111107
    112108void G4eCoulombScatteringModel::Initialise(const G4ParticleDefinition* p,
    113                                            const G4DataVector&)
    114 {
    115   //  G4cout << "!!! G4eCoulombScatteringModel::Initialise for "
    116   // << p->GetParticleName() << "  cos(TetMin)= " << cosThetaMin
    117   // << "  cos(TetMax)= " << cosThetaMax <<G4endl;
     109                                           const G4DataVector& cuts)
     110{
     111  SetupParticle(p);
     112  currentCouple = 0;
     113  elecXSection = nucXSection = 0.0;
     114  tkin = targetZ = mom2 = DBL_MIN;
     115  ecut = etag = DBL_MAX;
     116  cosThetaMin = cos(PolarAngleLimit());
     117  currentCuts = &cuts;
     118  //G4cout << "!!! G4eCoulombScatteringModel::Initialise for "
     119  //     << p->GetParticleName() << "  cos(TetMin)= " << cosThetaMin
     120  //     << "  cos(TetMax)= " << cosThetaMax <<G4endl;
    118121  if(!isInitialised) {
    119122    isInitialised = true;
    … …  
    124127    else
    125128      fParticleChange = new G4ParticleChangeForGamma();
    126   } else {
    127     return;
    128   }
    129 
    130   if(p->GetParticleType() == "nucleus") buildTable = false;
    131   if(!buildTable) return;
    132 
    133   // Compute log cross section table per atom
    134   if(!theCrossSectionTable) theCrossSectionTable = new G4PhysicsTable();
    135  
    136   nbins = 2*G4int(log10(highKEnergy/lowKEnergy));
    137 }
    138 
    139 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    140 
    141 G4double G4eCoulombScatteringModel::ComputeCrossSectionPerAtom(
    142                 const G4ParticleDefinition* p,
    143                 G4double kinEnergy,
    144                 G4double Z, G4double A,
    145                 G4double cutEnergy, G4double)
    146 {
    147   if(p == particle && kinEnergy == tkin && Z == targetZ &&
    148      A == targetA && cutEnergy == ecut) return nucXSection;
    149 
    150   //G4cout << "### G4eCoulombScatteringModel::ComputeCrossSectionPerAtom  for "
    151   //     << p->GetParticleName() << " Z= " << Z << " A= " << A
    152   //     << " e= " << kinEnergy << G4endl;
    153 
    154   nucXSection = ComputeElectronXSectionPerAtom(p,kinEnergy,Z,A,cutEnergy);
    155 
    156   // nuclear cross section
    157   if(theCrossSectionTable) {
    158     G4bool b;
    159     G4int iz  = G4int(Z);
    160     G4int idx = index[iz];
    161 
    162     // compute table for given Z
    163     if(-1 == idx) {
    164       idx = theCrossSectionTable->size();
    165       index[iz] = idx;
    166       G4PhysicsLogVector* ptrVector
    167         = new G4PhysicsLogVector(lowKEnergy, highKEnergy, nbins);
    168       //  G4cout << "New vector Z= " << iz << " A= " << A << " idx= " << idx << G4endl;
    169       G4double e, value;
    170       for(G4int i=0; i<=nbins; i++) {
    171         e     = ptrVector->GetLowEdgeEnergy( i ) ;
    172         value = CalculateCrossSectionPerAtom(p, e, Z, A); 
    173         ptrVector->PutValue( i, log(value) );
    174       }
    175       theCrossSectionTable->push_back(ptrVector);
    176     }
    177 
    178       // take value from the table
    179     nucXSection +=
    180       std::exp((((*theCrossSectionTable)[idx]))->GetValue(kinEnergy, b));
    181 
    182     // compute value from scratch
    183   } else nucXSection += CalculateCrossSectionPerAtom(p, kinEnergy, Z, A);
    184  
    185   //  G4cout << " cross(bn)= " << nucXSection/barn << G4endl;
    186  
    187   if(nucXSection < 0.0) nucXSection = 0.0;
    188   return nucXSection;
    189 }
    190 
    191 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    192 
    193 G4double G4eCoulombScatteringModel::ComputeElectronXSectionPerAtom(
    194                 const G4ParticleDefinition* p,
    195                 G4double kinEnergy,
    196                 G4double Z,
    197                 G4double A,
    198                 G4double cutEnergy)
    199 {
    200   if(p == particle && kinEnergy == tkin && Z == targetZ &&
    201      cutEnergy == ecut) return elecXSection;
     129  }
     130  if(mass < GeV && particle->GetParticleType() != "nucleus") {
     131    InitialiseElementSelectors(p,cuts);
     132  }
     133}
     134
     135//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     136
     137void G4eCoulombScatteringModel::ComputeMaxElectronScattering(G4double cutEnergy)
     138{
    202139  ecut = cutEnergy;
    203   elecXSection = 0.0;
    204   SetupParticle(p);
    205   G4double ekin = std::max(keV, kinEnergy);
    206   //G4double ekin = kinEnergy;
    207   SetupTarget(Z, A, ekin);
    208 
    209140  G4double tmax = tkin;
    210   if(p == theElectron) tmax *= 0.5;
    211   else if(p != thePositron) {
     141  cosTetMaxElec = 1.0;
     142  if(mass > MeV) {
    212143    G4double ratio = electron_mass_c2/mass;
    213144    G4double tau = tkin/mass;
    214145    tmax = 2.0*electron_mass_c2*tau*(tau + 2.)/
    215146      (1.0 + 2.0*ratio*(tau + 1.0) + ratio*ratio);
    216   }
    217 
    218   cosTetMaxElec = cosTetMaxNuc;
    219   G4double t = std::min(cutEnergy, tmax);
    220   G4double mom21 = t*(t + 2.0*electron_mass_c2);
    221   G4double t1 = tkin - t;
    222   if(t1 > 0.0) {
    223     G4double mom22 = t1*(t1 + 2.0*mass);
    224     G4double ctm = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22);
    225     if(ctm > cosTetMaxElec && ctm <= 1.0) cosTetMaxElec = ctm;
    226   }
    227 
    228   if(cosTetMaxElec < cosThetaMin) {
    229     G4double x1 = 1.0 - cosThetaMin  + screenZ;
    230     G4double x2 = 1.0 - cosTetMaxElec + screenZ;
    231     elecXSection = coeff*Z*chargeSquare*invbeta2*
    232       (cosThetaMin - cosTetMaxElec)/(x1*x2*mom2);
    233   }
    234   //  G4cout << "cut= " << ecut << " e= " << tkin
    235   // << " croosE(barn)= " << elecXSection/barn
    236   // << " cosEl= " << cosTetMaxElec << " costmin= " << cosThetaMin << G4endl;
    237   return elecXSection;
    238 }
    239 
    240 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    241 
    242 G4double G4eCoulombScatteringModel::CalculateCrossSectionPerAtom(
    243                              const G4ParticleDefinition* p,
    244                              G4double kinEnergy,
    245                              G4double Z, G4double A)
    246 {
    247   G4double cross = 0.0;
     147    cosTetMaxElec = 1.0 - std::min(cutEnergy, tmax)*electron_mass_c2/mom2;
     148  } else {
     149
     150    if(particle == theElectron) tmax *= 0.5;
     151    G4double t = std::min(cutEnergy, tmax);
     152    G4double mom21 = t*(t + 2.0*electron_mass_c2);
     153    G4double t1 = tkin - t;
     154    //G4cout << "tkin= " << tkin << " t= " << t << " t1= " << t1 << G4endl;
     155    if(t1 > 0.0) {
     156      G4double mom22 = t1*(t1 + 2.0*mass);
     157      G4double ctm = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22);
     158      //G4cout << "ctm= " << ctm << G4endl;
     159      if(ctm < 1.0) cosTetMaxElec = ctm;
     160    }
     161  }
     162}
     163
     164//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     165
     166G4double G4eCoulombScatteringModel::ComputeCrossSectionPerAtom(
     167                const G4ParticleDefinition* p,
     168                G4double kinEnergy,
     169                G4double Z, G4double,
     170                G4double cutEnergy, G4double)
     171{
     172  //G4cout << "### G4eCoulombScatteringModel::ComputeCrossSectionPerAtom  for "
     173  //  << p->GetParticleName()<<" Z= "<<Z<<" e(MeV)= "<< kinEnergy/MeV << G4endl;
     174  G4double xsec = 0.0;
    248175  SetupParticle(p);
    249   G4double ekin = std::max(keV, kinEnergy);
    250   //G4double ekin = kinEnergy;
    251   SetupTarget(Z, A, ekin);
    252 
    253   if(cosTetMaxNuc < cosThetaMin) {
    254     G4double x1 = 1.0 - cosThetaMin;
    255     G4double x2 = 1.0 - cosTetMaxNuc;
    256     G4double x3 = cosThetaMin - cosTetMaxNuc;
    257     G4double z1 = x1 + screenZ;
    258     G4double z2 = x2 + screenZ;
    259     G4double d  = 1.0/formfactA - screenZ;
    260     G4double d1 = 1.0 - formfactA*screenZ;
    261     G4double zn1= x1 + d;
    262     G4double zn2= x2 + d;
    263     cross = coeff*Z*Z*chargeSquare*invbeta2
    264       *(x3/(z1*z2) + x3/(zn1*zn2) +
    265         2.0*std::log(z1*zn2/(z2*zn1))/d) / (mom2*d1*d1);
    266   }
     176  G4double ekin = std::max(lowEnergyLimit, kinEnergy);
     177  SetupKinematic(ekin, cutEnergy);
     178  if(cosTetMaxNuc < cosTetMinNuc) {
     179    SetupTarget(Z, ekin);
     180    xsec = CrossSectionPerAtom(); 
     181  }
     182  /*
     183  G4cout << "e(MeV)= " << ekin/MeV << "cosTetMinNuc= " << cosTetMinNuc
     184         << " cosTetMaxNuc= " << cosTetMaxNuc
     185         << " cosTetMaxElec= " << cosTetMaxElec
     186         << " screenZ= " << screenZ
     187         << " formfactA= " << formfactA
     188         << " cosTetMaxHad= " << cosTetMaxHad << G4endl;
     189  */
     190  return xsec; 
     191}
     192
     193//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     194
     195G4double G4eCoulombScatteringModel::CrossSectionPerAtom()
     196{
     197  // This method needs initialisation before be called
     198
     199  G4double fac = coeff*targetZ*chargeSquare*invbeta2/mom2;
     200  elecXSection = 0.0;
     201  nucXSection  = 0.0;
     202
     203  G4double x  = 1.0 - cosTetMinNuc;
     204  G4double x1 = x + screenZ;
     205
     206  if(cosTetMaxElec2 < cosTetMinNuc) {
     207    elecXSection = fac*(cosTetMinNuc - cosTetMaxElec2)/
     208      (x1*(1.0 - cosTetMaxElec2 + screenZ));
     209    nucXSection  = elecXSection;
     210  }
     211
     212  //G4cout << "XS tkin(MeV)= " << tkin<<" xs= " <<nucXSection
     213  //     << " costmax= " << cosTetMaxNuc2
     214  //     << " costmin= " << cosTetMinNuc << "  Z= " << targetZ <<G4endl;
     215  if(cosTetMaxNuc2 < cosTetMinNuc) {
     216    G4double s  = screenZ*formfactA;
     217    G4double z1 = 1.0 - cosTetMaxNuc2 + screenZ;
     218    G4double d  = (1.0 - s)/formfactA;
     219    //G4cout <<"x1= "<<x1<<" z1= " <<z1<<" s= "<<s << " d= " <<d <<G4endl;
     220    if(d < 0.2*x1) {
     221      G4double x2 = x1*x1;
     222      G4double z2 = z1*z1;
     223      x = (1.0/(x1*x2) - 1.0/(z1*z2) - d*1.5*(1.0/(x2*x2) - 1.0/(z2*z2)))/
     224        (3.0*formfactA*formfactA);
     225    } else {
     226      G4double x2 = x1 + d;
     227      G4double z2 = z1 + d;
     228      x = (1.0 + 2.0*s)*((cosTetMinNuc - cosTetMaxNuc2)*(1.0/(x1*z1) + 1.0/(x2*z2)) -
     229                         2.0*log(z1*x2/(z2*x1))/d);
     230    }
     231    nucXSection += fac*targetZ*x;
     232  }
     233
     234  //G4cout<<" cross(bn)= "<<nucXSection/barn<<" xsElec(bn)= "<<elecXSection/barn
     235  //    << " Asc= " << screenZ << G4endl;
    267236 
    268   //  G4cout << "CalculateCrossSectionPerAtom: e(MeV)= " << tkin
    269   // << " cross(b)= " << cross/barn << " ctmin= " << cosThetaMin
    270   // << " ctmax= " << cosTetMaxNuc << G4endl;
    271  
    272   return cross;
     237  return nucXSection;
    273238}
    274239
    … …  
    276241
    277242void G4eCoulombScatteringModel::SampleSecondaries(
    278                 std::vector<G4DynamicParticle*>*,
     243                std::vector<G4DynamicParticle*>* fvect,
    279244                const G4MaterialCutsCouple* couple,
    280245                const G4DynamicParticle* dp,
    281246                G4double cutEnergy,
    282                 G4double maxEnergy)
    283 {
    284   const G4Material* aMaterial = couple->GetMaterial();
    285   const G4ParticleDefinition* p = dp->GetDefinition();
     247                G4double)
     248{
    286249  G4double kinEnergy = dp->GetKineticEnergy();
    287 
    288   // Select atom and setup
    289   SetupParticle(p);
    290   const G4Element* elm =
    291     SelectRandomAtom(aMaterial,p,kinEnergy,cutEnergy,maxEnergy);
    292   G4double Z  = elm->GetZ();
    293   G4double A  = elm->GetN();
    294 
    295   G4double cross =
    296     ComputeCrossSectionPerAtom(p,kinEnergy,Z,A,cutEnergy,maxEnergy);
    297 
    298   G4double costm = cosTetMaxNuc;
     250  if(kinEnergy <= DBL_MIN) return;
     251  DefineMaterial(couple);
     252  SetupParticle(dp->GetDefinition());
     253  G4double ekin = std::max(lowEnergyLimit, kinEnergy);
     254  SetupKinematic(ekin, cutEnergy);
     255  //G4cout << "G4eCoulombScatteringModel::SampleSecondaries e(MeV)= "
     256  //     << kinEnergy << "  " << particle->GetParticleName() << G4endl;
     257 
     258  // Choose nucleus
     259  currentElement = SelectRandomAtom(couple,particle,ekin,cutEnergy,ekin);
     260
     261  SetupTarget(currentElement->GetZ(),ekin);
     262 
     263  G4double cost = SampleCosineTheta();
     264  G4double z1   = 1.0 - cost;
     265  if(z1 < 0.0) return;
     266
     267  G4double sint = sqrt(z1*(1.0 + cost));
     268 
     269  //G4cout<<"## Sampled sint= " << sint << "  Z= " << targetZ
     270  //    << "  screenZ= " << screenZ << " cn= " << formfactA << G4endl;
     271 
     272  G4double phi  = twopi * G4UniformRand();
     273
     274  G4ThreeVector direction = dp->GetMomentumDirection();
     275  G4ThreeVector newDirection(cos(phi)*sint,sin(phi)*sint,cost);
     276  newDirection.rotateUz(direction);   
     277
     278  fParticleChange->ProposeMomentumDirection(newDirection);   
     279
     280  // recoil sampling assuming a small recoil
     281  // and first order correction to primary 4-momentum
     282  if(lowEnergyLimit < kinEnergy) {
     283    G4int ia = SelectIsotopeNumber(currentElement);
     284    G4double Trec = z1*mom2/(amu_c2*G4double(ia));
     285    G4double th =
     286      std::min(recoilThreshold,
     287               targetZ*currentElement->GetIonisation()->GetMeanExcitationEnergy());
     288
     289    if(Trec > th) {
     290      G4int iz = G4int(targetZ);
     291      G4ParticleDefinition* ion = theParticleTable->FindIon(iz, ia, 0, iz);
     292      Trec = z1*mom2/ion->GetPDGMass();
     293      if(Trec < kinEnergy) {
     294        G4ThreeVector dir = (direction - newDirection).unit();
     295        G4DynamicParticle* newdp = new G4DynamicParticle(ion, dir, Trec);
     296        fvect->push_back(newdp);
     297        fParticleChange->SetProposedKineticEnergy(kinEnergy - Trec);
     298      }
     299    }
     300  }
     301 
     302  return;
     303}
     304
     305//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     306
     307G4double G4eCoulombScatteringModel::SampleCosineTheta()
     308{
     309  G4double costm = cosTetMaxNuc2;
    299310  G4double formf = formfactA;
    300   if(G4UniformRand()*cross < elecXSection) {
    301     costm = cosTetMaxElec;
     311  G4double prob  = 0.0;
     312  G4double xs = CrossSectionPerAtom();
     313  if(xs > 0.0) prob = elecXSection/xs;
     314
     315  // scattering off e or A?
     316  if(G4UniformRand() < prob) {
     317    costm = cosTetMaxElec2;
    302318    formf = 0.0;
    303319  }
     320
    304321  /*
    305   G4cout << "G4eCoul...SampleSecondaries: e(MeV)= " << tkin
     322  G4cout << "SampleCost: e(MeV)= " << tkin
    306323         << " ctmin= " << cosThetaMin
    307324         << " ctmaxN= " << cosTetMaxNuc
    308325         << " ctmax= " << costm
    309          << " Z= " << Z << " A= " << A
    310          << " cross= " << cross/barn << " crossE= " << elecXSection/barn
     326         << " Z= " << targetZ << " A= " << targetA
    311327         << G4endl;
    312328  */
    313   if(costm >= cosThetaMin) return;
    314 
    315   G4double x1 = 1. - cosThetaMin + screenZ;
    316   G4double x2 = 1. - costm;
    317   G4double x3 = cosThetaMin - costm;
    318   G4double grej,  z, z1;
     329  if(costm >= cosTetMinNuc) return 2.0;
     330
     331  G4double x1 = 1. - cosTetMinNuc + screenZ;
     332  G4double x2 = 1. - costm + screenZ;
     333  G4double x3 = cosTetMinNuc - costm;
     334  G4double grej, z1;
    319335  do {
    320     z  = G4UniformRand()*x3;
    321     z1 = (x1*x2 - screenZ*z)/(x1 + z);
    322     if(z1 < 0.0) z1 = 0.0;
    323     else if(z1 > 2.0) z1 = 2.0;
     336    z1 = x1*x2/(x1 + G4UniformRand()*x3) - screenZ;
    324337    grej = 1.0/(1.0 + formf*z1);
    325338  } while ( G4UniformRand() > grej*grej ); 
    326  
    327   G4double cost = 1.0 - z1;
    328   G4double sint= sqrt(z1*(2.0 - z1));
    329   /*
    330   if(sint > 0.1)
    331     G4cout<<"## SampleSecondaries: e(MeV)= " << kinEnergy
    332           << " sint= " << sint << "  Z= " << Z << "  screenZ= " << screenZ
    333           << " cn= " << formf
    334           << G4endl;
    335   */
    336   G4double phi  = twopi * G4UniformRand();
    337 
    338   G4ThreeVector direction = dp->GetMomentumDirection();
    339   G4ThreeVector newDirection(cos(phi)*sint,sin(phi)*sint,cost);
    340   newDirection.rotateUz(direction);   
    341 
    342   fParticleChange->ProposeMomentumDirection(newDirection);   
    343  
    344   return;
    345 }
    346 
    347 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    348 
    349 
     339
     340  //G4cout << "z= " << z1 << " cross= " << nucXSection/barn
     341  // << " crossE= " << elecXSection/barn << G4endl;
     342
     343  return 1.0 - z1;
     344}
     345
     346//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     347
     348
Note: See TracChangeset for help on using the changeset viewer.