Geant4_10
G4AdjointeIonisationModel.cc
Go to the documentation of this file.
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 // $Id: G4AdjointeIonisationModel.cc 66892 2013-01-17 10:57:59Z gunter $
27 //
29 #include "G4AdjointCSManager.hh"
30 
31 #include "G4PhysicalConstants.hh"
32 #include "G4Integrator.hh"
33 #include "G4TrackStatus.hh"
34 #include "G4ParticleChange.hh"
35 #include "G4AdjointElectron.hh"
36 #include "G4Gamma.hh"
37 #include "G4AdjointGamma.hh"
38 
39 
41 //
43  G4VEmAdjointModel("Inv_eIon_model")
44 
45 {
46 
47  UseMatrix =true;
48  UseMatrixPerElement = true;
49  ApplyCutInRange = true;
52  WithRapidSampling = false;
53 
58 }
60 //
62 {;}
64 //
66  G4bool IsScatProjToProjCase,
67  G4ParticleChange* fParticleChange)
68 {
69 
70 
71  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
72 
73  //Elastic inverse scattering
74  //---------------------------------------------------------
75  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
76  G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
77 
78  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
79  return;
80  }
81 
82  //Sample secondary energy
83  //-----------------------
84  G4double projectileKinEnergy;
85  if (!WithRapidSampling ) { //used by default
86  projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, IsScatProjToProjCase);
87 
88  CorrectPostStepWeight(fParticleChange,
89  aTrack.GetWeight(),
90  adjointPrimKinEnergy,
91  projectileKinEnergy,
92  IsScatProjToProjCase); //Caution !!!this weight correction should be always applied
93  }
94  else { //only for test at the moment
95 
96  G4double Emin,Emax;
97  if (IsScatProjToProjCase) {
99  Emax=GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
100  }
101  else {
102  Emin=GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);
103  Emax=GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
104  }
105  projectileKinEnergy = Emin*std::pow(Emax/Emin,G4UniformRand());
106 
107 
108 
110  if ( !IsScatProjToProjCase) lastCS=lastAdjointCSForProdToProjCase;
111 
112  G4double new_weight=aTrack.GetWeight();
113  G4double used_diffCS=lastCS*std::log(Emax/Emin)/projectileKinEnergy;
114  G4double needed_diffCS=adjointPrimKinEnergy/projectileKinEnergy;
115  if (!IsScatProjToProjCase) needed_diffCS *=DiffCrossSectionPerVolumePrimToSecond(currentMaterial,projectileKinEnergy,adjointPrimKinEnergy);
116  else needed_diffCS *=DiffCrossSectionPerVolumePrimToScatPrim(currentMaterial,projectileKinEnergy,adjointPrimKinEnergy);
117  new_weight*=needed_diffCS/used_diffCS;
118  fParticleChange->SetParentWeightByProcess(false);
119  fParticleChange->SetSecondaryWeightByProcess(false);
120  fParticleChange->ProposeParentWeight(new_weight);
121 
122 
123  }
124 
125 
126 
127  //Kinematic:
128  //we consider a two body elastic scattering for the forward processes where the projectile knock on an e- at rest and gives
129  // him part of its energy
130  //----------------------------------------------------------------------------------------
131 
133  G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
134  G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
135 
136 
137 
138  //Companion
139  //-----------
141  if (IsScatProjToProjCase) {
143  }
144  G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
145  G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;
146 
147 
148  //Projectile momentum
149  //--------------------
150  G4double P_parallel = (adjointPrimP*adjointPrimP + projectileP2 - companionP2)/(2.*adjointPrimP);
151  G4double P_perp = std::sqrt( projectileP2 - P_parallel*P_parallel);
152  G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
153  G4double phi =G4UniformRand()*2.*3.1415926;
154  G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
155  projectileMomentum.rotateUz(dir_parallel);
156 
157 
158 
159  if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
160  fParticleChange->ProposeTrackStatus(fStopAndKill);
161  fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
162  //G4cout<<"projectileMomentum "<<projectileMomentum<<G4endl;
163  }
164  else {
165  fParticleChange->ProposeEnergy(projectileKinEnergy);
166  fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
167  }
168 
169 
170 
171 
172 }
174 //
175 //The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
177  G4double kinEnergyProj,
178  G4double kinEnergyProd,
179  G4double Z,
180  G4double )
181 {
182  G4double dSigmadEprod=0;
183  G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
184  G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
185 
186 
187  if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ //the produced particle should have a kinetic energy smaller than the projectile
188  dSigmadEprod=Z*DiffCrossSectionMoller(kinEnergyProj,kinEnergyProd);
189  }
190  return dSigmadEprod;
191 
192 
193 
194 }
195 
197 //
198 G4double G4AdjointeIonisationModel::DiffCrossSectionMoller(G4double kinEnergyProj,G4double kinEnergyProd){
199 
200  G4double energy = kinEnergyProj + electron_mass_c2;
201  G4double x = kinEnergyProd/kinEnergyProj;
202  G4double gam = energy/electron_mass_c2;
203  G4double gamma2 = gam*gam;
204  G4double beta2 = 1.0 - 1.0/gamma2;
205 
206  G4double gg = (2.0*gam - 1.0)/gamma2;
207  G4double y = 1.0 - x;
209  G4double dCS = fac*( 1.-gg + ((1.0 - gg*x)/(x*x))
210  + ((1.0 - gg*y)/(y*y)))/(beta2*(gam-1));
211  return dCS/kinEnergyProj;
212 
213 
214 
215 }
216 
virtual void SampleSecondaries(const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)
virtual G4double DiffCrossSectionPerVolumePrimToSecond(const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyProd)
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase(G4double PrimAdjEnergy)
G4double lastAdjointCSForScatProjToProjCase
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
const G4DynamicParticle * GetDynamicParticle() const
G4double SampleAdjSecEnergyFromCSMatrix(size_t MatrixIndex, G4double prim_energy, G4bool IsScatProjToProjCase)
G4ParticleDefinition * theDirectPrimaryPartDef
void ProposeParentWeight(G4double finalWeight)
virtual G4double DiffCrossSectionPerVolumePrimToScatPrim(const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyScatProj)
static G4AdjointElectron * AdjointElectron()
tuple x
Definition: test.py:50
G4Material * currentMaterial
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
G4double GetTotalMomentum() const
Double_t y
Definition: plot.C:279
virtual G4double DiffCrossSectionPerAtomPrimToSecond(G4double kinEnergyProj, G4double kinEnergyProd, G4double Z, G4double A=0.)
void SetSecondaryWeightByProcess(G4bool)
void SetParentWeightByProcess(G4bool)
#define G4UniformRand()
Definition: Randomize.hh:87
double energy
Definition: plottest35.C:25
virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
Float_t Z
Definition: plot.C:39
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:72
float electron_mass_c2
Definition: hepunit.py:274
G4bool UseOnlyOneMatrixForAllElements
G4double GetPDGMass() const
virtual G4double GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy, G4double Tcut=0)
Hep3Vector unit() const
void ProposeEnergy(G4double finalEnergy)
void AddSecondary(G4Track *aSecondary)
G4double GetWeight() const
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4ParticleDefinition * theAdjEquivOfDirectSecondPartDef
G4double currentTcutForDirectSecond
G4double lastAdjointCSForProdToProjCase
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
virtual void CorrectPostStepWeight(G4ParticleChange *fParticleChange, G4double old_weight, G4double adjointPrimKinEnergy, G4double projectileKinEnergy, G4bool IsScatProjToProjCase)