| 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 #include "G4AdjointCSManager.hh"
28 29
29 #include "G4AdjointElectron.hh" <<
30 #include "G4Electron.hh" <<
31 #include "G4ParticleChange.hh" <<
32 #include "G4PhysicalConstants.hh" 30 #include "G4PhysicalConstants.hh"
>> 31 #include "G4Integrator.hh"
33 #include "G4TrackStatus.hh" 32 #include "G4TrackStatus.hh"
>> 33 #include "G4ParticleChange.hh"
>> 34 #include "G4AdjointElectron.hh"
>> 35 #include "G4Gamma.hh"
>> 36 #include "G4AdjointGamma.hh"
>> 37
34 38
35 ////////////////////////////////////////////// 39 ////////////////////////////////////////////////////////////////////////////////
36 G4AdjointeIonisationModel::G4AdjointeIonisatio << 40 //
37 : G4VEmAdjointModel("Inv_eIon_model") << 41 G4AdjointeIonisationModel::G4AdjointeIonisationModel():
>> 42 G4VEmAdjointModel("Inv_eIon_model")
38 43
39 { 44 {
40 fUseMatrix = true; << 45
41 fUseMatrixPerElement = true; << 46 UseMatrix =true;
42 fApplyCutInRange = true; << 47 UseMatrixPerElement = true;
43 fOneMatrixForAllElements = true; << 48 ApplyCutInRange = true;
44 << 49 UseOnlyOneMatrixForAllElements = true;
45 fAdjEquivDirectPrimPart = G4AdjointElectro << 50 CS_biasing_factor =1.;
46 fAdjEquivDirectSecondPart = G4AdjointElectro << 51 WithRapidSampling = false;
47 fDirectPrimaryPart = G4Electron::Elec << 52
48 fSecondPartSameType = true; << 53 theAdjEquivOfDirectPrimPartDef =G4AdjointElectron::AdjointElectron();
>> 54 theAdjEquivOfDirectSecondPartDef=G4AdjointElectron::AdjointElectron();
>> 55 theDirectPrimaryPartDef=G4Electron::Electron();
>> 56 second_part_of_same_type=true;
49 } 57 }
50 <<
51 ////////////////////////////////////////////// 58 ////////////////////////////////////////////////////////////////////////////////
52 G4AdjointeIonisationModel::~G4AdjointeIonisati << 59 //
53 << 60 G4AdjointeIonisationModel::~G4AdjointeIonisationModel()
>> 61 {;}
54 ////////////////////////////////////////////// 62 ////////////////////////////////////////////////////////////////////////////////
55 void G4AdjointeIonisationModel::SampleSecondar << 63 //
56 const G4Track& aTrack, G4bool IsScatProjToPr << 64 void G4AdjointeIonisationModel::SampleSecondaries(const G4Track& aTrack,
57 G4ParticleChange* fParticleChange) << 65 G4bool IsScatProjToProjCase,
58 { << 66 G4ParticleChange* fParticleChange)
59 const G4DynamicParticle* theAdjointPrimary = << 67 {
>> 68
>> 69
>> 70 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
>> 71
>> 72 //Elastic inverse scattering
>> 73 //---------------------------------------------------------
>> 74 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
>> 75 G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
>> 76
>> 77 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
>> 78 return;
>> 79 }
>> 80
>> 81 //Sample secondary energy
>> 82 //-----------------------
>> 83 G4double projectileKinEnergy;
>> 84 if (!WithRapidSampling ) { //used by default
>> 85 projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, IsScatProjToProjCase);
>> 86
>> 87 CorrectPostStepWeight(fParticleChange,
>> 88 aTrack.GetWeight(),
>> 89 adjointPrimKinEnergy,
>> 90 projectileKinEnergy,
>> 91 IsScatProjToProjCase); //Caution !!!this weight correction should be always applied
>> 92 }
>> 93 else { //only for test at the moment
>> 94
>> 95 G4double Emin,Emax;
>> 96 if (IsScatProjToProjCase) {
>> 97 Emin=GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond);
>> 98 Emax=GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
>> 99 }
>> 100 else {
>> 101 Emin=GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);
>> 102 Emax=GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
>> 103 }
>> 104 projectileKinEnergy = Emin*std::pow(Emax/Emin,G4UniformRand());
>> 105
>> 106
>> 107
>> 108 lastCS=lastAdjointCSForScatProjToProjCase;
>> 109 if ( !IsScatProjToProjCase) lastCS=lastAdjointCSForProdToProjCase;
>> 110
>> 111 G4double new_weight=aTrack.GetWeight();
>> 112 G4double used_diffCS=lastCS*std::log(Emax/Emin)/projectileKinEnergy;
>> 113 G4double needed_diffCS=adjointPrimKinEnergy/projectileKinEnergy;
>> 114 if (!IsScatProjToProjCase) needed_diffCS *=DiffCrossSectionPerVolumePrimToSecond(currentMaterial,projectileKinEnergy,adjointPrimKinEnergy);
>> 115 else needed_diffCS *=DiffCrossSectionPerVolumePrimToScatPrim(currentMaterial,projectileKinEnergy,adjointPrimKinEnergy);
>> 116 new_weight*=needed_diffCS/used_diffCS;
>> 117 fParticleChange->SetParentWeightByProcess(false);
>> 118 fParticleChange->SetSecondaryWeightByProcess(false);
>> 119 fParticleChange->ProposeParentWeight(new_weight);
>> 120
>> 121
>> 122 }
>> 123
>> 124
>> 125
>> 126 //Kinematic:
>> 127 //we consider a two body elastic scattering for the forward processes where the projectile knock on an e- at rest and gives
>> 128 // him part of its energy
>> 129 //----------------------------------------------------------------------------------------
>> 130
>> 131 G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
>> 132 G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
>> 133 G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
>> 134
>> 135
>> 136
>> 137 //Companion
>> 138 //-----------
>> 139 G4double companionM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
>> 140 if (IsScatProjToProjCase) {
>> 141 companionM0=theAdjEquivOfDirectSecondPartDef->GetPDGMass();
>> 142 }
>> 143 G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
>> 144 G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;
>> 145
>> 146
>> 147 //Projectile momentum
>> 148 //--------------------
>> 149 G4double P_parallel = (adjointPrimP*adjointPrimP + projectileP2 - companionP2)/(2.*adjointPrimP);
>> 150 G4double P_perp = std::sqrt( projectileP2 - P_parallel*P_parallel);
>> 151 G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
>> 152 G4double phi =G4UniformRand()*2.*3.1415926;
>> 153 G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
>> 154 projectileMomentum.rotateUz(dir_parallel);
>> 155
>> 156
>> 157
>> 158 if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
>> 159 fParticleChange->ProposeTrackStatus(fStopAndKill);
>> 160 fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
>> 161 //G4cout<<"projectileMomentum "<<projectileMomentum<<G4endl;
>> 162 }
>> 163 else {
>> 164 fParticleChange->ProposeEnergy(projectileKinEnergy);
>> 165 fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
>> 166 }
>> 167
60 168
61 // Elastic inverse scattering << 169
62 G4double adjointPrimKinEnergy = theAdjointPr <<
63 G4double adjointPrimP = theAdjointPr <<
64 <<
65 if(adjointPrimKinEnergy > GetHighEnergyLimit <<
66 { <<
67 return; <<
68 } <<
69 <<
70 // Sample secondary energy <<
71 G4double projectileKinEnergy; <<
72 if(!fWithRapidSampling) <<
73 { // used by default <<
74 projectileKinEnergy = <<
75 SampleAdjSecEnergyFromCSMatrix(adjointPr <<
76 <<
77 // Caution!!! this weight correction shoul <<
78 CorrectPostStepWeight(fParticleChange, aTr <<
79 adjointPrimKinEnergy <<
80 IsScatProjToProj); <<
81 } <<
82 else <<
83 { // only for testing <<
84 G4double Emin, Emax; <<
85 if(IsScatProjToProj) <<
86 { <<
87 Emin = GetSecondAdjEnergyMinForScatProjT <<
88 <<
89 Emax = GetSecondAdjEnergyMaxForScatProjT <<
90 } <<
91 else <<
92 { <<
93 Emin = GetSecondAdjEnergyMinForProdToPro <<
94 Emax = GetSecondAdjEnergyMaxForProdToPro <<
95 } <<
96 projectileKinEnergy = Emin * std::pow(Emax <<
97 <<
98 fLastCS = fLastAdjointCSForScatProjToProj; <<
99 if(!IsScatProjToProj) <<
100 fLastCS = fLastAdjointCSForProdToProj; <<
101 <<
102 G4double new_weight = aTrack.GetWeight(); <<
103 G4double used_diffCS = <<
104 fLastCS * std::log(Emax / Emin) / projec <<
105 G4double needed_diffCS = adjointPrimKinEne <<
106 if(!IsScatProjToProj) <<
107 needed_diffCS *= DiffCrossSectionPerVolu <<
108 fCurrentMaterial, projectileKinEnergy, <<
109 else <<
110 needed_diffCS *= DiffCrossSectionPerVolu <<
111 fCurrentMaterial, projectileKinEnergy, <<
112 new_weight *= needed_diffCS / used_diffCS; <<
113 fParticleChange->SetParentWeightByProcess( <<
114 fParticleChange->SetSecondaryWeightByProce <<
115 fParticleChange->ProposeParentWeight(new_w <<
116 } <<
117 <<
118 // Kinematic: <<
119 // we consider a two body elastic scattering <<
120 // the projectile knock on an e- at rest and <<
121 G4double projectileM0 = fAdjEquivDi <<
122 G4double projectileTotalEnergy = projectileM <<
123 G4double projectileP2 = <<
124 projectileTotalEnergy * projectileTotalEne <<
125 <<
126 // Companion <<
127 G4double companionM0 = fAdjEquivDirectPrimPa <<
128 if(IsScatProjToProj) <<
129 { <<
130 companionM0 = fAdjEquivDirectSecondPart->G <<
131 } <<
132 G4double companionTotalEnergy = <<
133 companionM0 + projectileKinEnergy - adjoin <<
134 G4double companionP2 = <<
135 companionTotalEnergy * companionTotalEnerg <<
136 <<
137 // Projectile momentum <<
138 G4double P_parallel = <<
139 (adjointPrimP * adjointPrimP + projectileP <<
140 (2. * adjointPrimP); <<
141 G4double P_perp = std::sqrt(projectileP2 - P <<
142 G4ThreeVector dir_parallel = theAdjointPrima <<
143 G4double phi = G4UniformRand() <<
144 G4ThreeVector projectileMomentum = <<
145 G4ThreeVector(P_perp * std::cos(phi), P_pe <<
146 projectileMomentum.rotateUz(dir_parallel); <<
147 <<
148 if(!IsScatProjToProj) <<
149 { // kill the primary and add a secondary <<
150 fParticleChange->ProposeTrackStatus(fStopA <<
151 fParticleChange->AddSecondary( <<
152 new G4DynamicParticle(fAdjEquivDirectPri <<
153 } <<
154 else <<
155 { <<
156 fParticleChange->ProposeEnergy(projectileK <<
157 fParticleChange->ProposeMomentumDirection( <<
158 } <<
159 } <<
160 170
>> 171 }
161 ////////////////////////////////////////////// 172 ////////////////////////////////////////////////////////////////////////////////
162 // The implementation here is correct for ener << 173 //
163 // photoelectric and compton scattering the me << 174 //The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
164 G4double G4AdjointeIonisationModel::DiffCrossS 175 G4double G4AdjointeIonisationModel::DiffCrossSectionPerAtomPrimToSecond(
165 G4double kinEnergyProj, G4double kinEnergyPr << 176 G4double kinEnergyProj,
>> 177 G4double kinEnergyProd,
>> 178 G4double Z,
>> 179 G4double )
166 { 180 {
167 G4double dSigmadEprod = 0.; << 181 G4double dSigmadEprod=0;
168 G4double Emax_proj = GetSecondAdjEnergyMa << 182 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
169 G4double Emin_proj = GetSecondAdjEnergyMi << 183 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
170 << 184
171 // the produced particle should have a kinet << 185
172 // projectile << 186 if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ //the produced particle should have a kinetic energy smaller than the projectile
173 if(kinEnergyProj > Emin_proj && kinEnergyPro << 187 dSigmadEprod=Z*DiffCrossSectionMoller(kinEnergyProj,kinEnergyProd);
174 { << 188 }
175 dSigmadEprod = Z * DiffCrossSectionMoller( << 189 return dSigmadEprod;
176 } << 190
177 return dSigmadEprod; << 191
>> 192
178 } 193 }
179 194
180 ////////////////////////////////////////////// 195 //////////////////////////////////////////////////////////////////////////////
181 G4double G4AdjointeIonisationModel::DiffCrossS << 196 //
182 G4double kinEnergyProj, G4double kinEnergyPr << 197 G4double G4AdjointeIonisationModel::DiffCrossSectionMoller(G4double kinEnergyProj,G4double kinEnergyProd){
183 { << 198
184 // G4double energy = kinEnergyProj + electro << 199 G4double energy = kinEnergyProj + electron_mass_c2;
185 G4double x = kinEnergyProd / kinEnergyP << 200 G4double x = kinEnergyProd/kinEnergyProj;
186 G4double gam = (kinEnergyProj + electron_ << 201 G4double gam = energy/electron_mass_c2;
187 G4double gamma2 = gam * gam; << 202 G4double gamma2 = gam*gam;
188 G4double beta2 = 1.0 - 1.0 / gamma2; << 203 G4double beta2 = 1.0 - 1.0/gamma2;
189 << 204
190 G4double gg = (2.0 * gam - 1.0) / gamma2; << 205 G4double gg = (2.0*gam - 1.0)/gamma2;
191 G4double y = 1.0 - x; << 206 G4double y = 1.0 - x;
192 G4double fac = twopi_mc2_rcl2 / electron_mas << 207 G4double fac=twopi_mc2_rcl2/electron_mass_c2;
193 G4double dCS = << 208 G4double dCS = fac*( 1.-gg + ((1.0 - gg*x)/(x*x))
194 fac * (1. - gg + ((1.0 - gg * x) / (x * x) << 209 + ((1.0 - gg*y)/(y*y)))/(beta2*(gam-1));
195 (beta2 * (gam - 1.)); << 210 return dCS/kinEnergyProj;
196 return dCS / kinEnergyProj; << 211
197 } << 212
>> 213
>> 214 }
>> 215
198 216