| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 26
27 #include "G4AdjointeIonisationModel.hh" 27 #include "G4AdjointeIonisationModel.hh"
28 28
29 #include "G4AdjointElectron.hh" 29 #include "G4AdjointElectron.hh"
30 #include "G4Electron.hh" 30 #include "G4Electron.hh"
31 #include "G4ParticleChange.hh" 31 #include "G4ParticleChange.hh"
32 #include "G4PhysicalConstants.hh" 32 #include "G4PhysicalConstants.hh"
33 #include "G4TrackStatus.hh" 33 #include "G4TrackStatus.hh"
34 34
35 ////////////////////////////////////////////// 35 ////////////////////////////////////////////////////////////////////////////////
36 G4AdjointeIonisationModel::G4AdjointeIonisatio 36 G4AdjointeIonisationModel::G4AdjointeIonisationModel()
37 : G4VEmAdjointModel("Inv_eIon_model") 37 : G4VEmAdjointModel("Inv_eIon_model")
38 38
39 { 39 {
40 fUseMatrix = true; 40 fUseMatrix = true;
41 fUseMatrixPerElement = true; 41 fUseMatrixPerElement = true;
42 fApplyCutInRange = true; 42 fApplyCutInRange = true;
43 fOneMatrixForAllElements = true; 43 fOneMatrixForAllElements = true;
44 44
45 fAdjEquivDirectPrimPart = G4AdjointElectro 45 fAdjEquivDirectPrimPart = G4AdjointElectron::AdjointElectron();
46 fAdjEquivDirectSecondPart = G4AdjointElectro 46 fAdjEquivDirectSecondPart = G4AdjointElectron::AdjointElectron();
47 fDirectPrimaryPart = G4Electron::Elec 47 fDirectPrimaryPart = G4Electron::Electron();
48 fSecondPartSameType = true; 48 fSecondPartSameType = true;
49 } 49 }
50 50
51 ////////////////////////////////////////////// 51 ////////////////////////////////////////////////////////////////////////////////
52 G4AdjointeIonisationModel::~G4AdjointeIonisati 52 G4AdjointeIonisationModel::~G4AdjointeIonisationModel() {}
53 53
54 ////////////////////////////////////////////// 54 ////////////////////////////////////////////////////////////////////////////////
55 void G4AdjointeIonisationModel::SampleSecondar 55 void G4AdjointeIonisationModel::SampleSecondaries(
56 const G4Track& aTrack, G4bool IsScatProjToPr 56 const G4Track& aTrack, G4bool IsScatProjToProj,
57 G4ParticleChange* fParticleChange) 57 G4ParticleChange* fParticleChange)
58 { 58 {
59 const G4DynamicParticle* theAdjointPrimary = 59 const G4DynamicParticle* theAdjointPrimary = aTrack.GetDynamicParticle();
60 60
61 // Elastic inverse scattering 61 // Elastic inverse scattering
62 G4double adjointPrimKinEnergy = theAdjointPr 62 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
63 G4double adjointPrimP = theAdjointPr 63 G4double adjointPrimP = theAdjointPrimary->GetTotalMomentum();
64 64
65 if(adjointPrimKinEnergy > GetHighEnergyLimit 65 if(adjointPrimKinEnergy > GetHighEnergyLimit() * 0.999)
66 { 66 {
67 return; 67 return;
68 } 68 }
69 69
70 // Sample secondary energy 70 // Sample secondary energy
71 G4double projectileKinEnergy; 71 G4double projectileKinEnergy;
72 if(!fWithRapidSampling) 72 if(!fWithRapidSampling)
73 { // used by default 73 { // used by default
74 projectileKinEnergy = 74 projectileKinEnergy =
75 SampleAdjSecEnergyFromCSMatrix(adjointPr 75 SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, IsScatProjToProj);
76 76
77 // Caution!!! this weight correction shoul 77 // Caution!!! this weight correction should be always applied
78 CorrectPostStepWeight(fParticleChange, aTr 78 CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(),
79 adjointPrimKinEnergy 79 adjointPrimKinEnergy, projectileKinEnergy,
80 IsScatProjToProj); 80 IsScatProjToProj);
81 } 81 }
82 else 82 else
83 { // only for testing 83 { // only for testing
84 G4double Emin, Emax; 84 G4double Emin, Emax;
85 if(IsScatProjToProj) 85 if(IsScatProjToProj)
86 { 86 {
87 Emin = GetSecondAdjEnergyMinForScatProjT 87 Emin = GetSecondAdjEnergyMinForScatProjToProj(adjointPrimKinEnergy,
88 88 fTcutSecond);
89 Emax = GetSecondAdjEnergyMaxForScatProjT 89 Emax = GetSecondAdjEnergyMaxForScatProjToProj(adjointPrimKinEnergy);
90 } 90 }
91 else 91 else
92 { 92 {
93 Emin = GetSecondAdjEnergyMinForProdToPro 93 Emin = GetSecondAdjEnergyMinForProdToProj(adjointPrimKinEnergy);
94 Emax = GetSecondAdjEnergyMaxForProdToPro 94 Emax = GetSecondAdjEnergyMaxForProdToProj(adjointPrimKinEnergy);
95 } 95 }
96 projectileKinEnergy = Emin * std::pow(Emax 96 projectileKinEnergy = Emin * std::pow(Emax / Emin, G4UniformRand());
97 97
98 fLastCS = fLastAdjointCSForScatProjToProj; 98 fLastCS = fLastAdjointCSForScatProjToProj;
99 if(!IsScatProjToProj) 99 if(!IsScatProjToProj)
100 fLastCS = fLastAdjointCSForProdToProj; 100 fLastCS = fLastAdjointCSForProdToProj;
101 101
102 G4double new_weight = aTrack.GetWeight(); 102 G4double new_weight = aTrack.GetWeight();
103 G4double used_diffCS = 103 G4double used_diffCS =
104 fLastCS * std::log(Emax / Emin) / projec 104 fLastCS * std::log(Emax / Emin) / projectileKinEnergy;
105 G4double needed_diffCS = adjointPrimKinEne 105 G4double needed_diffCS = adjointPrimKinEnergy / projectileKinEnergy;
106 if(!IsScatProjToProj) 106 if(!IsScatProjToProj)
107 needed_diffCS *= DiffCrossSectionPerVolu 107 needed_diffCS *= DiffCrossSectionPerVolumePrimToSecond(
108 fCurrentMaterial, projectileKinEnergy, 108 fCurrentMaterial, projectileKinEnergy, adjointPrimKinEnergy);
109 else 109 else
110 needed_diffCS *= DiffCrossSectionPerVolu 110 needed_diffCS *= DiffCrossSectionPerVolumePrimToScatPrim(
111 fCurrentMaterial, projectileKinEnergy, 111 fCurrentMaterial, projectileKinEnergy, adjointPrimKinEnergy);
112 new_weight *= needed_diffCS / used_diffCS; 112 new_weight *= needed_diffCS / used_diffCS;
113 fParticleChange->SetParentWeightByProcess( 113 fParticleChange->SetParentWeightByProcess(false);
114 fParticleChange->SetSecondaryWeightByProce 114 fParticleChange->SetSecondaryWeightByProcess(false);
115 fParticleChange->ProposeParentWeight(new_w 115 fParticleChange->ProposeParentWeight(new_weight);
116 } 116 }
117 117
118 // Kinematic: 118 // Kinematic:
119 // we consider a two body elastic scattering 119 // we consider a two body elastic scattering for the forward processes where
120 // the projectile knock on an e- at rest and 120 // the projectile knock on an e- at rest and gives it part of its energy
121 G4double projectileM0 = fAdjEquivDi 121 G4double projectileM0 = fAdjEquivDirectPrimPart->GetPDGMass();
122 G4double projectileTotalEnergy = projectileM 122 G4double projectileTotalEnergy = projectileM0 + projectileKinEnergy;
123 G4double projectileP2 = 123 G4double projectileP2 =
124 projectileTotalEnergy * projectileTotalEne 124 projectileTotalEnergy * projectileTotalEnergy - projectileM0 * projectileM0;
125 125
126 // Companion 126 // Companion
127 G4double companionM0 = fAdjEquivDirectPrimPa 127 G4double companionM0 = fAdjEquivDirectPrimPart->GetPDGMass();
128 if(IsScatProjToProj) 128 if(IsScatProjToProj)
129 { 129 {
130 companionM0 = fAdjEquivDirectSecondPart->G 130 companionM0 = fAdjEquivDirectSecondPart->GetPDGMass();
131 } 131 }
132 G4double companionTotalEnergy = 132 G4double companionTotalEnergy =
133 companionM0 + projectileKinEnergy - adjoin 133 companionM0 + projectileKinEnergy - adjointPrimKinEnergy;
134 G4double companionP2 = 134 G4double companionP2 =
135 companionTotalEnergy * companionTotalEnerg 135 companionTotalEnergy * companionTotalEnergy - companionM0 * companionM0;
136 136
137 // Projectile momentum 137 // Projectile momentum
138 G4double P_parallel = 138 G4double P_parallel =
139 (adjointPrimP * adjointPrimP + projectileP 139 (adjointPrimP * adjointPrimP + projectileP2 - companionP2) /
140 (2. * adjointPrimP); 140 (2. * adjointPrimP);
141 G4double P_perp = std::sqrt(projectileP2 - P 141 G4double P_perp = std::sqrt(projectileP2 - P_parallel * P_parallel);
142 G4ThreeVector dir_parallel = theAdjointPrima 142 G4ThreeVector dir_parallel = theAdjointPrimary->GetMomentumDirection();
143 G4double phi = G4UniformRand() 143 G4double phi = G4UniformRand() * twopi;
144 G4ThreeVector projectileMomentum = 144 G4ThreeVector projectileMomentum =
145 G4ThreeVector(P_perp * std::cos(phi), P_pe 145 G4ThreeVector(P_perp * std::cos(phi), P_perp * std::sin(phi), P_parallel);
146 projectileMomentum.rotateUz(dir_parallel); 146 projectileMomentum.rotateUz(dir_parallel);
147 147
148 if(!IsScatProjToProj) 148 if(!IsScatProjToProj)
149 { // kill the primary and add a secondary 149 { // kill the primary and add a secondary
150 fParticleChange->ProposeTrackStatus(fStopA 150 fParticleChange->ProposeTrackStatus(fStopAndKill);
151 fParticleChange->AddSecondary( 151 fParticleChange->AddSecondary(
152 new G4DynamicParticle(fAdjEquivDirectPri 152 new G4DynamicParticle(fAdjEquivDirectPrimPart, projectileMomentum));
153 } 153 }
154 else 154 else
155 { 155 {
156 fParticleChange->ProposeEnergy(projectileK 156 fParticleChange->ProposeEnergy(projectileKinEnergy);
157 fParticleChange->ProposeMomentumDirection( 157 fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
158 } 158 }
159 } 159 }
160 160
161 ////////////////////////////////////////////// 161 ////////////////////////////////////////////////////////////////////////////////
162 // The implementation here is correct for ener 162 // The implementation here is correct for energy loss process, for the
163 // photoelectric and compton scattering the me 163 // photoelectric and compton scattering the method should be redefined
164 G4double G4AdjointeIonisationModel::DiffCrossS 164 G4double G4AdjointeIonisationModel::DiffCrossSectionPerAtomPrimToSecond(
165 G4double kinEnergyProj, G4double kinEnergyPr 165 G4double kinEnergyProj, G4double kinEnergyProd, G4double Z, G4double)
166 { 166 {
167 G4double dSigmadEprod = 0.; 167 G4double dSigmadEprod = 0.;
168 G4double Emax_proj = GetSecondAdjEnergyMa 168 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProj(kinEnergyProd);
169 G4double Emin_proj = GetSecondAdjEnergyMi 169 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProj(kinEnergyProd);
170 170
171 // the produced particle should have a kinet 171 // the produced particle should have a kinetic energy smaller than the
172 // projectile 172 // projectile
173 if(kinEnergyProj > Emin_proj && kinEnergyPro 173 if(kinEnergyProj > Emin_proj && kinEnergyProj <= Emax_proj)
174 { 174 {
175 dSigmadEprod = Z * DiffCrossSectionMoller( 175 dSigmadEprod = Z * DiffCrossSectionMoller(kinEnergyProj, kinEnergyProd);
176 } 176 }
177 return dSigmadEprod; 177 return dSigmadEprod;
178 } 178 }
179 179
180 ////////////////////////////////////////////// 180 //////////////////////////////////////////////////////////////////////////////
181 G4double G4AdjointeIonisationModel::DiffCrossS 181 G4double G4AdjointeIonisationModel::DiffCrossSectionMoller(
182 G4double kinEnergyProj, G4double kinEnergyPr 182 G4double kinEnergyProj, G4double kinEnergyProd)
183 { 183 {
184 // G4double energy = kinEnergyProj + electro 184 // G4double energy = kinEnergyProj + electron_mass_c2;
185 G4double x = kinEnergyProd / kinEnergyP 185 G4double x = kinEnergyProd / kinEnergyProj;
186 G4double gam = (kinEnergyProj + electron_ 186 G4double gam = (kinEnergyProj + electron_mass_c2) / electron_mass_c2;
187 G4double gamma2 = gam * gam; 187 G4double gamma2 = gam * gam;
188 G4double beta2 = 1.0 - 1.0 / gamma2; 188 G4double beta2 = 1.0 - 1.0 / gamma2;
189 189
190 G4double gg = (2.0 * gam - 1.0) / gamma2; 190 G4double gg = (2.0 * gam - 1.0) / gamma2;
191 G4double y = 1.0 - x; 191 G4double y = 1.0 - x;
192 G4double fac = twopi_mc2_rcl2 / electron_mas 192 G4double fac = twopi_mc2_rcl2 / electron_mass_c2;
193 G4double dCS = 193 G4double dCS =
194 fac * (1. - gg + ((1.0 - gg * x) / (x * x) 194 fac * (1. - gg + ((1.0 - gg * x) / (x * x)) + ((1.0 - gg * y) / (y * y))) /
195 (beta2 * (gam - 1.)); 195 (beta2 * (gam - 1.));
196 return dCS / kinEnergyProj; 196 return dCS / kinEnergyProj;
197 } 197 }
198 198