source: trunk/source/processes/electromagnetic/standard/src/G4IonCoulombScatteringModel.cc @ 1350

Last change on this file since 1350 was 1350, checked in by garnier, 14 years ago

update to last version 4.9.4
File size: 8.1 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer                                           *
4// *                                                                  *
5// * The  Geant4 software  is  copyright of the Copyright Holders  of *
6// * the Geant4 Collaboration.  It is provided  under  the terms  and *
7// * conditions of the Geant4 Software License,  included in the file *
8// * LICENSE and available at  http://cern.ch/geant4/license .  These *
9// * include a list of copyright holders.                             *
10// *                                                                  *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work  make  any representation or  warranty, express or implied, *
14// * regarding  this  software system or assume any liability for its *
15// * use.  Please see the license in the file  LICENSE  and URL above *
16// * for the full disclaimer and the limitation of liability.         *
17// *                                                                  *
18// * This  code  implementation is the result of  the  scientific and *
19// * technical work of the GEANT4 collaboration.                      *
20// * By using,  copying,  modifying or  distributing the software (or *
21// * any work based  on the software)  you  agree  to acknowledge its *
22// * use  in  resulting  scientific  publications,  and indicate your *
23// * acceptance of all terms of the Geant4 Software license.          *
24// ********************************************************************
25//
26//      G4IonCoulombScatteringModel.cc
27// -------------------------------------------------------------------
28//
29// GEANT4 Class header file
30//
31// File name:    G4IonCoulombScatteringModel
32//
33// Author:      Cristina Consolandi
34//
35// Creation date: 05.10.2010 from G4eCoulombScatteringModel
36//                               & G4CoulombScatteringModel
37//
38// Class Description:
39//      Single Scattering Model for
40//      for protons, alpha and heavy Ions
41//
42// Reference:
43//      M.J. Boschini et al. "Nuclear and Non-Ionizing Energy-Loss
44//      for Coulomb ScatteredParticles from Low Energy up to Relativistic
45//      Regime in Space Radiation Environment"
46//      Accepted for publication in the Proceedings of  the  ICATPP Conference
47//      on Cosmic Rays for Particle and Astroparticle Physics, Villa  Olmo, 7-8
48//      October,  2010, to be published by World Scientific (Singapore).
49//
50//      Available for downloading at:
51//      http://arxiv.org/abs/1011.4822
52//
53// -------------------------------------------------------------------
54//
55//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
56
57
58#include "G4IonCoulombScatteringModel.hh"
59#include "Randomize.hh"
60//#include "G4DataVector.hh"
61#include "G4ParticleChangeForGamma.hh"
62#include "G4Proton.hh"
63#include "G4ProductionCutsTable.hh"
64#include "G4NucleiProperties.hh"
65
66#include "G4UnitsTable.hh"
67
68//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
69
70using namespace std;
71
72G4IonCoulombScatteringModel::G4IonCoulombScatteringModel(const G4String& nam)
73  : G4VEmModel(nam),
74
75    cosThetaMin(1.0),
76    isInitialised(false)
77{
78        fNistManager = G4NistManager::Instance();
79        theParticleTable = G4ParticleTable::GetParticleTable();
80        theProton   = G4Proton::Proton();
81
82        pCuts=0;
83        currentMaterial = 0;
84        currentElement  = 0;
85        currentCouple = 0;
86
87        lowEnergyLimit  = 100*eV;
88        recoilThreshold = 0.*eV;
89        heavycorr =0;
90        particle = 0;
91        mass=0;
92
93        ioncross = new G4IonCoulombCrossSection(); 
94
95}
96
97
98//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
99
100G4IonCoulombScatteringModel::~G4IonCoulombScatteringModel()
101{ delete  ioncross;}
102
103//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
104
105void G4IonCoulombScatteringModel::Initialise(const G4ParticleDefinition* p,
106                                           const G4DataVector&  )
107{
108        SetupParticle(p);
109        currentCouple = 0;
110        cosThetaMin = cos(PolarAngleLimit());
111        ioncross->Initialise(p,cosThetaMin);
112 
113        pCuts = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(3);
114
115
116        if(!isInitialised) {
117                isInitialised = true;
118                fParticleChange = GetParticleChangeForGamma();
119                }
120}
121
122//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
123
124G4double G4IonCoulombScatteringModel::ComputeCrossSectionPerAtom(
125                                const G4ParticleDefinition* p,
126                                G4double kinEnergy, 
127                                G4double Z, 
128                                G4double, 
129                                G4double cutEnergy,
130                                G4double)
131{
132
133        SetupParticle(p);
134 
135        G4double xsec =0.0;
136                if(kinEnergy < lowEnergyLimit) return xsec;
137
138        DefineMaterial(CurrentCouple());
139
140        G4int iz          = G4int(Z);
141
142        //from lab to pCM & mu_rel of effective particle
143        ioncross->SetupKinematic(kinEnergy, cutEnergy,iz);
144
145
146        ioncross->SetupTarget(Z, kinEnergy, heavycorr);
147
148
149        xsec = ioncross->NuclearCrossSection();
150
151//cout<< "..........xsec "<<G4BestUnit(xsec,"Surface") <<endl;
152  return xsec;
153}
154
155//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
156
157void G4IonCoulombScatteringModel::SampleSecondaries(
158                               std::vector<G4DynamicParticle*>* fvect,
159                               const G4MaterialCutsCouple* couple,
160                               const G4DynamicParticle* dp,
161                               G4double cutEnergy, 
162                               G4double)
163{
164        G4double kinEnergy = dp->GetKineticEnergy();
165
166        if(kinEnergy < lowEnergyLimit) return;
167       
168        DefineMaterial(couple);
169
170        SetupParticle(dp->GetDefinition());
171
172        // Choose nucleus
173        currentElement = SelectRandomAtom(couple,particle,
174                                    kinEnergy,cutEnergy,kinEnergy);
175
176        G4double Z  = currentElement->GetZ();
177        G4int iz          = G4int(Z);
178        G4int ia = SelectIsotopeNumber(currentElement);
179        G4double m2 = G4NucleiProperties::GetNuclearMass(ia, iz);
180
181
182
183        G4double xsec= ComputeCrossSectionPerAtom(particle,kinEnergy, Z,
184                                kinEnergy, cutEnergy, kinEnergy) ;
185        if(xsec == 0.0)return;
186   
187        //scattering angle, z1 == (1-cost)
188        G4double z1 = ioncross->SampleCosineTheta(); 
189
190        if(z1 <= 0.0) { return; }
191
192        G4double cost = 1.0 - z1;
193        G4double sint = sqrt(z1*(1.0 + cost));
194        G4double phi  = twopi * G4UniformRand();
195
196
197        // kinematics in the Lab system
198        G4double etot = kinEnergy + mass;
199        G4double mom2=  kinEnergy*(kinEnergy+2.0*mass);
200        G4double ptot = sqrt(mom2);
201
202        //CM particle 1
203        G4double bet  = ptot/(etot + m2);
204        G4double gam  = 1.0/sqrt((1.0 - bet)*(1.0 + bet));
205
206        //CM   
207        G4double momCM2= ioncross->GetMomentum2();
208        G4double momCM =std::sqrt(momCM2); 
209        //energy & momentum after scattering of incident particle
210        G4double pxCM = momCM*sint*cos(phi);
211        G4double pyCM = momCM*sint*sin(phi);
212        G4double pzCM = momCM*cost;
213        G4double eCM  = sqrt(momCM2 + mass*mass);
214
215        //CM--->Lab
216        G4ThreeVector v1(pxCM , pyCM, gam*(pzCM + bet*eCM));
217        G4ThreeVector dir = dp->GetMomentumDirection(); 
218
219        G4ThreeVector newDirection = v1.unit();
220        newDirection.rotateUz(dir);   
221 
222        fParticleChange->ProposeMomentumDirection(newDirection);   
223 
224        // recoil.......................................
225        G4double trec =(1.0 - cost)* m2*(etot*etot - mass*mass )/
226                                (mass*mass + m2*m2+ 2.*m2*etot);
227
228        G4double finalT = kinEnergy - trec;
229
230
231        if(finalT <= lowEnergyLimit) { 
232                trec = kinEnergy; 
233                finalT = 0.0;
234                } 
235   
236        fParticleChange->SetProposedKineticEnergy(finalT);
237
238        G4double tcut = recoilThreshold;
239        if(pCuts) { tcut= std::max(tcut,(*pCuts)[currentMaterialIndex]); 
240
241                        //G4cout<<" tcut eV "<<tcut/eV<<endl;
242                        }
243 
244        if(trec > tcut) {
245                G4ParticleDefinition* ion = theParticleTable->FindIon(iz, ia, 0, iz);
246                G4double plab = sqrt(finalT*(finalT + 2.0*mass));
247                G4ThreeVector p2 = (ptot*dir - plab*newDirection).unit();
248                G4DynamicParticle* newdp  = new G4DynamicParticle(ion, p2, trec);
249                        fvect->push_back(newdp);
250        } else if(trec > 0.0) {
251                fParticleChange->ProposeLocalEnergyDeposit(trec);
252                fParticleChange->ProposeNonIonizingEnergyDeposit(trec);
253        }
254
255
256}
257
258//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
259               
Note: See TracBrowser for help on using the repository browser.