| 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 // $Id: G4AdjointeIonisationModel.cc 66892 2013-01-17 10:57:59Z gunter $
>> 27 //
27 #include "G4AdjointeIonisationModel.hh" 28 #include "G4AdjointeIonisationModel.hh"
>> 29 #include "G4AdjointCSManager.hh"
28 30
29 #include "G4AdjointElectron.hh" <<
30 #include "G4Electron.hh" <<
31 #include "G4ParticleChange.hh" <<
32 #include "G4PhysicalConstants.hh" 31 #include "G4PhysicalConstants.hh"
>> 32 #include "G4Integrator.hh"
33 #include "G4TrackStatus.hh" 33 #include "G4TrackStatus.hh"
>> 34 #include "G4ParticleChange.hh"
>> 35 #include "G4AdjointElectron.hh"
>> 36 #include "G4Gamma.hh"
>> 37 #include "G4AdjointGamma.hh"
>> 38
34 39
35 ////////////////////////////////////////////// 40 ////////////////////////////////////////////////////////////////////////////////
36 G4AdjointeIonisationModel::G4AdjointeIonisatio << 41 //
37 : G4VEmAdjointModel("Inv_eIon_model") << 42 G4AdjointeIonisationModel::G4AdjointeIonisationModel():
>> 43 G4VEmAdjointModel("Inv_eIon_model")
38 44
39 { 45 {
40 fUseMatrix = true; << 46
41 fUseMatrixPerElement = true; << 47 UseMatrix =true;
42 fApplyCutInRange = true; << 48 UseMatrixPerElement = true;
43 fOneMatrixForAllElements = true; << 49 ApplyCutInRange = true;
44 << 50 UseOnlyOneMatrixForAllElements = true;
45 fAdjEquivDirectPrimPart = G4AdjointElectro << 51 CS_biasing_factor =1.;
46 fAdjEquivDirectSecondPart = G4AdjointElectro << 52 WithRapidSampling = false;
47 fDirectPrimaryPart = G4Electron::Elec << 53
48 fSecondPartSameType = true; << 54 theAdjEquivOfDirectPrimPartDef =G4AdjointElectron::AdjointElectron();
>> 55 theAdjEquivOfDirectSecondPartDef=G4AdjointElectron::AdjointElectron();
>> 56 theDirectPrimaryPartDef=G4Electron::Electron();
>> 57 second_part_of_same_type=true;
49 } 58 }
50 <<
51 ////////////////////////////////////////////// 59 ////////////////////////////////////////////////////////////////////////////////
52 G4AdjointeIonisationModel::~G4AdjointeIonisati << 60 //
53 << 61 G4AdjointeIonisationModel::~G4AdjointeIonisationModel()
>> 62 {;}
54 ////////////////////////////////////////////// 63 ////////////////////////////////////////////////////////////////////////////////
55 void G4AdjointeIonisationModel::SampleSecondar << 64 //
56 const G4Track& aTrack, G4bool IsScatProjToPr << 65 void G4AdjointeIonisationModel::SampleSecondaries(const G4Track& aTrack,
57 G4ParticleChange* fParticleChange) << 66 G4bool IsScatProjToProjCase,
58 { << 67 G4ParticleChange* fParticleChange)
59 const G4DynamicParticle* theAdjointPrimary = << 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) {
>> 98 Emin=GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond);
>> 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
>> 109 lastCS=lastAdjointCSForScatProjToProjCase;
>> 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
>> 132 G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
>> 133 G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
>> 134 G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
>> 135
>> 136
>> 137
>> 138 //Companion
>> 139 //-----------
>> 140 G4double companionM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
>> 141 if (IsScatProjToProjCase) {
>> 142 companionM0=theAdjEquivOfDirectSecondPartDef->GetPDGMass();
>> 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
60 169
61 // Elastic inverse scattering << 170
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 171
>> 172 }
161 ////////////////////////////////////////////// 173 ////////////////////////////////////////////////////////////////////////////////
162 // The implementation here is correct for ener << 174 //
163 // photoelectric and compton scattering the me << 175 //The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
164 G4double G4AdjointeIonisationModel::DiffCrossS 176 G4double G4AdjointeIonisationModel::DiffCrossSectionPerAtomPrimToSecond(
165 G4double kinEnergyProj, G4double kinEnergyPr << 177 G4double kinEnergyProj,
>> 178 G4double kinEnergyProd,
>> 179 G4double Z,
>> 180 G4double )
166 { 181 {
167 G4double dSigmadEprod = 0.; << 182 G4double dSigmadEprod=0;
168 G4double Emax_proj = GetSecondAdjEnergyMa << 183 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
169 G4double Emin_proj = GetSecondAdjEnergyMi << 184 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
170 << 185
171 // the produced particle should have a kinet << 186
172 // projectile << 187 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 << 188 dSigmadEprod=Z*DiffCrossSectionMoller(kinEnergyProj,kinEnergyProd);
174 { << 189 }
175 dSigmadEprod = Z * DiffCrossSectionMoller( << 190 return dSigmadEprod;
176 } << 191
177 return dSigmadEprod; << 192
>> 193
178 } 194 }
179 195
180 ////////////////////////////////////////////// 196 //////////////////////////////////////////////////////////////////////////////
181 G4double G4AdjointeIonisationModel::DiffCrossS << 197 //
182 G4double kinEnergyProj, G4double kinEnergyPr << 198 G4double G4AdjointeIonisationModel::DiffCrossSectionMoller(G4double kinEnergyProj,G4double kinEnergyProd){
183 { << 199
184 // G4double energy = kinEnergyProj + electro << 200 G4double energy = kinEnergyProj + electron_mass_c2;
185 G4double x = kinEnergyProd / kinEnergyP << 201 G4double x = kinEnergyProd/kinEnergyProj;
186 G4double gam = (kinEnergyProj + electron_ << 202 G4double gam = energy/electron_mass_c2;
187 G4double gamma2 = gam * gam; << 203 G4double gamma2 = gam*gam;
188 G4double beta2 = 1.0 - 1.0 / gamma2; << 204 G4double beta2 = 1.0 - 1.0/gamma2;
189 << 205
190 G4double gg = (2.0 * gam - 1.0) / gamma2; << 206 G4double gg = (2.0*gam - 1.0)/gamma2;
191 G4double y = 1.0 - x; << 207 G4double y = 1.0 - x;
192 G4double fac = twopi_mc2_rcl2 / electron_mas << 208 G4double fac=twopi_mc2_rcl2/electron_mass_c2;
193 G4double dCS = << 209 G4double dCS = fac*( 1.-gg + ((1.0 - gg*x)/(x*x))
194 fac * (1. - gg + ((1.0 - gg * x) / (x * x) << 210 + ((1.0 - gg*y)/(y*y)))/(beta2*(gam-1));
195 (beta2 * (gam - 1.)); << 211 return dCS/kinEnergyProj;
196 return dCS / kinEnergyProj; << 212
197 } << 213
>> 214
>> 215 }
>> 216
198 217