G4eCoulombScatteringModel.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028 // -------------------------------------------------------------------
00029 //
00030 // GEANT4 Class file
00031 //
00032 //
00033 // File name:     G4eCoulombScatteringModel
00034 //
00035 // Author:        Vladimir Ivanchenko 
00036 //
00037 // Creation date: 22.08.2005
00038 //
00039 // Modifications:
00040 //
00041 // 01.08.06 V.Ivanchenko extend upper limit of table to TeV and review the
00042 //          logic of building - only elements from G4ElementTable
00043 // 08.08.06 V.Ivanchenko build internal table in ekin scale, introduce faclim
00044 // 19.08.06 V.Ivanchenko add inline function ScreeningParameter 
00045 // 09.10.07 V.Ivanchenko reorganized methods, add cut dependence in scattering off e- 
00046 // 09.06.08 V.Ivanchenko add SelectIsotope and sampling of the recoil ion 
00047 // 16.06.09 C.Consolandi fixed computation of effective mass
00048 // 27.05.10 V.Ivanchenko added G4WentzelOKandVIxSection class to
00049 //              compute cross sections and sample scattering angle
00050 //
00051 //
00052 // Class Description:
00053 //
00054 // -------------------------------------------------------------------
00055 //
00056 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00057 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00058 
00059 #include "G4eCoulombScatteringModel.hh"
00060 #include "G4PhysicalConstants.hh"
00061 #include "G4SystemOfUnits.hh"
00062 #include "Randomize.hh"
00063 #include "G4DataVector.hh"
00064 #include "G4ElementTable.hh"
00065 #include "G4ParticleChangeForGamma.hh"
00066 #include "G4Proton.hh"
00067 #include "G4ParticleTable.hh"
00068 #include "G4ProductionCutsTable.hh"
00069 #include "G4NucleiProperties.hh"
00070 #include "G4Pow.hh"
00071 #include "G4LossTableManager.hh"
00072 #include "G4LossTableBuilder.hh"
00073 #include "G4NistManager.hh"
00074 
00075 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00076 
00077 using namespace std;
00078 
00079 G4eCoulombScatteringModel::G4eCoulombScatteringModel(const G4String& nam)
00080   : G4VEmModel(nam),
00081     cosThetaMin(1.0),
00082     cosThetaMax(-1.0),
00083     isInitialised(false)
00084 {
00085   fParticleChange = 0;
00086   fNistManager = G4NistManager::Instance();
00087   theParticleTable = G4ParticleTable::GetParticleTable();
00088   theProton   = G4Proton::Proton();
00089   currentMaterial = 0; 
00090 
00091   pCuts = 0;
00092 
00093   lowEnergyThreshold = 1*keV;  // particle will be killed for lower energy
00094   recoilThreshold = 0.*keV; // by default does not work
00095 
00096   particle = 0;
00097   currentCouple = 0;
00098   wokvi = new G4WentzelOKandVIxSection();
00099 
00100   currentMaterialIndex = 0;
00101 
00102   cosTetMinNuc = 1.0;
00103   cosTetMaxNuc = -1.0;
00104   elecRatio = 0.0;
00105   mass = proton_mass_c2;
00106 }
00107 
00108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00109 
00110 G4eCoulombScatteringModel::~G4eCoulombScatteringModel()
00111 {
00112   delete wokvi;
00113 }
00114 
00115 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00116 
00117 void G4eCoulombScatteringModel::Initialise(const G4ParticleDefinition* p,
00118                                            const G4DataVector& cuts)
00119 {
00120   SetupParticle(p);
00121   currentCouple = 0;
00122   cosThetaMin = cos(PolarAngleLimit());
00123   wokvi->Initialise(p, cosThetaMin);
00124   /*    
00125   G4cout << "G4eCoulombScatteringModel: " << particle->GetParticleName()
00126          << "  1-cos(ThetaLimit)= " << 1 - cosThetaMin
00127          << "  cos(thetaMax)= " <<  cosThetaMax
00128          << G4endl;
00129   */
00130   pCuts = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(3);
00131   /*
00132   G4cout << "!!! G4eCoulombScatteringModel::Initialise for " 
00133          << p->GetParticleName() << "  cos(TetMin)= " << cosThetaMin 
00134          << "  cos(TetMax)= " << cosThetaMax <<G4endl;
00135   G4cout << "cut0= " << cuts[0] << "  cut1= " << cuts[1] << G4endl;
00136   */
00137   if(!isInitialised) {
00138     isInitialised = true;
00139     fParticleChange = GetParticleChangeForGamma();
00140   }
00141   if(mass < GeV) {
00142     InitialiseElementSelectors(p,cuts);
00143   }
00144 }
00145 
00146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00147 
00148 G4double G4eCoulombScatteringModel::ComputeCrossSectionPerAtom(
00149                 const G4ParticleDefinition* p,
00150                 G4double kinEnergy,
00151                 G4double Z, G4double,
00152                 G4double cutEnergy, G4double)
00153 {
00154   //G4cout << "### G4eCoulombScatteringModel::ComputeCrossSectionPerAtom  for " 
00155   //  << p->GetParticleName()<<" Z= "<<Z<<" e(MeV)= "<< kinEnergy/MeV << G4endl; 
00156   G4double cross = 0.0;
00157   if(p != particle) { SetupParticle(p); }
00158 
00159   // cross section is set to zero to avoid problems in sample secondary
00160   if(kinEnergy <= 0.0) { return cross; }
00161   DefineMaterial(CurrentCouple());
00162   cosTetMinNuc = wokvi->SetupKinematic(kinEnergy, currentMaterial);
00163   if(cosThetaMax < cosTetMinNuc) {
00164     G4int iz = G4int(Z);
00165     cosTetMinNuc = wokvi->SetupTarget(iz, cutEnergy);
00166     cosTetMaxNuc = cosThetaMax; 
00167     if(iz == 1 && cosTetMaxNuc < 0.0 && particle == theProton) { 
00168       cosTetMaxNuc = 0.0; 
00169     }
00170     cross =  wokvi->ComputeNuclearCrossSection(cosTetMinNuc, cosTetMaxNuc);
00171     elecRatio = wokvi->ComputeElectronCrossSection(cosTetMinNuc, cosThetaMax);
00172     cross += elecRatio;
00173     if(cross > 0.0) { elecRatio /= cross; }  
00174   }
00175   /*
00176   if(p->GetParticleName() == "e-") 
00177   G4cout << "e(MeV)= " << kinEnergy/MeV << " cross(b)= " << cross/barn  
00178          << " 1-cosTetMinNuc= " << 1-cosTetMinNuc
00179          << " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc
00180          << G4endl;
00181   */
00182   return cross;  
00183 }
00184 
00185 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00186 
00187 void G4eCoulombScatteringModel::SampleSecondaries(
00188                 std::vector<G4DynamicParticle*>* fvect,
00189                 const G4MaterialCutsCouple* couple,
00190                 const G4DynamicParticle* dp,
00191                 G4double cutEnergy,
00192                 G4double)
00193 {
00194   G4double kinEnergy = dp->GetKineticEnergy();
00195 
00196   // absorb particle below low-energy limit to avoid situation
00197   // when a particle has no energy loss
00198   if(kinEnergy < lowEnergyThreshold) { 
00199     fParticleChange->SetProposedKineticEnergy(0.0);
00200     fParticleChange->ProposeLocalEnergyDeposit(kinEnergy);
00201     fParticleChange->ProposeNonIonizingEnergyDeposit(kinEnergy);
00202     return; 
00203   }
00204   SetupParticle(dp->GetDefinition());
00205   DefineMaterial(couple);
00206   /*
00207   G4cout << "G4eCoulombScatteringModel::SampleSecondaries e(MeV)= " 
00208          << kinEnergy << "  " << particle->GetParticleName() 
00209          << " cut= " << cutEnergy<< G4endl;
00210   */
00211   // Choose nucleus
00212   const G4Element* currentElement = 
00213     SelectRandomAtom(couple,particle,kinEnergy,cutEnergy,kinEnergy);
00214 
00215   G4double Z = currentElement->GetZ();
00216 
00217   if(ComputeCrossSectionPerAtom(particle,kinEnergy, Z,
00218                                 kinEnergy, cutEnergy, kinEnergy) == 0.0) 
00219     { return; }
00220 
00221   G4int iz = G4int(Z);
00222   G4int ia = SelectIsotopeNumber(currentElement);
00223   G4double targetMass = G4NucleiProperties::GetNuclearMass(ia, iz);
00224   wokvi->SetTargetMass(targetMass);
00225 
00226   G4ThreeVector newDirection = 
00227     wokvi->SampleSingleScattering(cosTetMinNuc, cosThetaMax, elecRatio);
00228   G4double cost = newDirection.z();
00229 
00230   G4ThreeVector direction = dp->GetMomentumDirection(); 
00231   newDirection.rotateUz(direction);   
00232 
00233   fParticleChange->ProposeMomentumDirection(newDirection);   
00234 
00235   // recoil sampling assuming a small recoil
00236   // and first order correction to primary 4-momentum
00237   G4double mom2 = wokvi->GetMomentumSquare();
00238   G4double trec = mom2*(1.0 - cost)/(targetMass + (mass + kinEnergy)*(1.0 - cost));
00239   G4double finalT = kinEnergy - trec; 
00240   //G4cout<<"G4eCoulombScatteringModel: finalT= "<<finalT<<" Trec= "<<trec<<G4endl;
00241   if(finalT <= lowEnergyThreshold) { 
00242     trec = kinEnergy;  
00243     finalT = 0.0;
00244   } 
00245 
00246   fParticleChange->SetProposedKineticEnergy(finalT);
00247   G4double tcut = recoilThreshold;
00248   if(pCuts) { tcut= std::max(tcut,(*pCuts)[currentMaterialIndex]); }
00249 
00250   if(trec > tcut) {
00251     G4ParticleDefinition* ion = theParticleTable->GetIon(iz, ia, 0.0);
00252     G4ThreeVector dir = (direction*sqrt(mom2) - 
00253                          newDirection*sqrt(finalT*(2*mass + finalT))).unit();
00254     G4DynamicParticle* newdp = new G4DynamicParticle(ion, dir, trec);
00255     fvect->push_back(newdp);
00256   } else {
00257     fParticleChange->ProposeLocalEnergyDeposit(trec);
00258     fParticleChange->ProposeNonIonizingEnergyDeposit(trec);
00259   }
00260  
00261   return;
00262 }
00263 
00264 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00265 
00266 

Generated on Mon May 27 17:48:05 2013 for Geant4 by  doxygen 1.4.7