| 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: G4AdjointIonIonisationModel.cc,v 1.3 2010-11-11 11:51:56 ldesorgh Exp $
>> 27 // GEANT4 tag $Name: not supported by cvs2svn $
>> 28 //
27 #include "G4AdjointIonIonisationModel.hh" 29 #include "G4AdjointIonIonisationModel.hh"
28 <<
29 #include "G4AdjointCSManager.hh" 30 #include "G4AdjointCSManager.hh"
>> 31
>> 32
>> 33 #include "G4Integrator.hh"
>> 34 #include "G4TrackStatus.hh"
>> 35 #include "G4ParticleChange.hh"
30 #include "G4AdjointElectron.hh" 36 #include "G4AdjointElectron.hh"
>> 37 #include "G4AdjointProton.hh"
>> 38 #include "G4AdjointInterpolator.hh"
31 #include "G4BetheBlochModel.hh" 39 #include "G4BetheBlochModel.hh"
32 #include "G4BraggIonModel.hh" 40 #include "G4BraggIonModel.hh"
>> 41 #include "G4Proton.hh"
33 #include "G4GenericIon.hh" 42 #include "G4GenericIon.hh"
34 #include "G4NistManager.hh" 43 #include "G4NistManager.hh"
35 #include "G4ParticleChange.hh" <<
36 #include "G4PhysicalConstants.hh" <<
37 #include "G4SystemOfUnits.hh" <<
38 #include "G4TrackStatus.hh" <<
39 44
40 ////////////////////////////////////////////// 45 ////////////////////////////////////////////////////////////////////////////////
41 G4AdjointIonIonisationModel::G4AdjointIonIonis << 46 //
42 : G4VEmAdjointModel("Adjoint_IonIonisation") << 47 G4AdjointIonIonisationModel::G4AdjointIonIonisationModel():
43 { << 48 G4VEmAdjointModel("Adjoint_IonIonisation")
44 fUseMatrix = true; << 49 {
45 fUseMatrixPerElement = true; << 50
46 fApplyCutInRange = true; << 51
47 fOneMatrixForAllElements = true; << 52 UseMatrix =true;
48 fSecondPartSameType = false; << 53 UseMatrixPerElement = true;
49 << 54 ApplyCutInRange = true;
50 // The direct EM Model is taken as BetheBloc << 55 UseOnlyOneMatrixForAllElements = true;
51 // computation of the differential cross sec << 56 CS_biasing_factor =1.;
52 // The Bragg model could be used as an alter << 57 second_part_of_same_type =false;
53 // differential cross section << 58 use_only_bragg = false; // for the Ion ionisation using the parametrised table model the cross sections and the sample of secondaries is done
54 << 59 // as in the BraggIonModel, Therefore the use of this flag;
55 fBetheBlochDirectEMModel = new G4BetheBloch << 60
56 fBraggIonDirectEMModel = new G4BraggIonMo << 61 //The direct EM Model is taken has BetheBloch it is only used for the computation
57 fAdjEquivDirectSecondPart = G4AdjointElectro << 62 // of the differential cross section.
58 fDirectPrimaryPart = nullptr; << 63 //The Bragg model could be used as an alternative as it offers the same differential cross section
59 } << 64
>> 65 theBetheBlochDirectEMModel = new G4BetheBlochModel(G4GenericIon::GenericIon());
>> 66 theBraggIonDirectEMModel = new G4BraggIonModel(G4GenericIon::GenericIon());
>> 67 theAdjEquivOfDirectSecondPartDef=G4AdjointElectron::AdjointElectron();
>> 68 theDirectPrimaryPartDef =0;
>> 69 theAdjEquivOfDirectPrimPartDef =0;
>> 70 /* theDirectPrimaryPartDef =fwd_ion;
>> 71 theAdjEquivOfDirectPrimPartDef =adj_ion;
>> 72
>> 73 DefineProjectileProperty();*/
60 74
61 ////////////////////////////////////////////// <<
62 G4AdjointIonIonisationModel::~G4AdjointIonIoni <<
63 75
64 ////////////////////////////////////////////// <<
65 void G4AdjointIonIonisationModel::SampleSecond <<
66 const G4Track& aTrack, G4bool isScatProjToPr <<
67 G4ParticleChange* fParticleChange) <<
68 { <<
69 const G4DynamicParticle* theAdjointPrimary = <<
70 76
71 // Elastic inverse scattering <<
72 G4double adjointPrimKinEnergy = theAdjointPr <<
73 G4double adjointPrimP = theAdjointPr <<
74 <<
75 if(adjointPrimKinEnergy > GetHighEnergyLimit <<
76 { <<
77 return; <<
78 } <<
79 77
80 // Sample secondary energy <<
81 G4double projectileKinEnergy = <<
82 SampleAdjSecEnergyFromCSMatrix(adjointPrim <<
83 // Caution !!!this weight correction should <<
84 CorrectPostStepWeight(fParticleChange, aTrac <<
85 adjointPrimKinEnergy, <<
86 isScatProjToProj); <<
87 <<
88 // Kinematics: <<
89 // we consider a two body elastic scattering <<
90 // the projectile knock on an e- at rest and <<
91 G4double projectileM0 = fAdjEquivDi <<
92 G4double projectileTotalEnergy = projectileM <<
93 G4double projectileP2 = <<
94 projectileTotalEnergy * projectileTotalEne <<
95 <<
96 // Companion <<
97 G4double companionM0 = fAdjEquivDirectPrimPa <<
98 if(isScatProjToProj) <<
99 { <<
100 companionM0 = fAdjEquivDirectSecondPart->G <<
101 } <<
102 G4double companionTotalEnergy = <<
103 companionM0 + projectileKinEnergy - adjoin <<
104 G4double companionP2 = <<
105 companionTotalEnergy * companionTotalEnerg <<
106 <<
107 // Projectile momentum <<
108 G4double P_parallel = <<
109 (adjointPrimP * adjointPrimP + projectileP <<
110 (2. * adjointPrimP); <<
111 G4double P_perp = std::sqrt(projectileP2 - P <<
112 G4ThreeVector dir_parallel = theAdjointPrima <<
113 G4double phi = G4UniformRand() <<
114 G4ThreeVector projectileMomentum = <<
115 G4ThreeVector(P_perp * std::cos(phi), P_pe <<
116 projectileMomentum.rotateUz(dir_parallel); <<
117 <<
118 if(!isScatProjToProj) <<
119 { // kill the primary and add a secondary <<
120 fParticleChange->ProposeTrackStatus(fStopA <<
121 fParticleChange->AddSecondary( <<
122 new G4DynamicParticle(fAdjEquivDirectPri <<
123 } <<
124 else <<
125 { <<
126 fParticleChange->ProposeEnergy(projectileK <<
127 fParticleChange->ProposeMomentumDirection( <<
128 } <<
129 } 78 }
130 <<
131 ////////////////////////////////////////////// 79 ////////////////////////////////////////////////////////////////////////////////
132 G4double G4AdjointIonIonisationModel::DiffCros << 80 //
133 G4double kinEnergyProj, G4double kinEnergyPr << 81 G4AdjointIonIonisationModel::~G4AdjointIonIonisationModel()
134 { << 82 {;}
135 // Probably that here the Bragg Model should << 83 ////////////////////////////////////////////////////////////////////////////////
136 // kinEnergyProj/nuc<2MeV << 84 //
137 G4double dSigmadEprod = 0.; << 85 void G4AdjointIonIonisationModel::SampleSecondaries(const G4Track& aTrack,
138 G4double Emax_proj = GetSecondAdjEnergyMa << 86 G4bool IsScatProjToProjCase,
139 G4double Emin_proj = GetSecondAdjEnergyMi << 87 G4ParticleChange* fParticleChange)
140 << 88 {
141 G4double kinEnergyProjScaled = fMassRatio * << 89 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
142 << 90
143 // the produced particle should have a kinet << 91 //Elastic inverse scattering
144 // projectile << 92 //---------------------------------------------------------
145 if(kinEnergyProj > Emin_proj && kinEnergyPro << 93 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
146 { << 94 G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
147 G4double Tmax = kinEnergyProj; << 95
148 << 96 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
149 G4double E1 = kinEnergyProd; << 97 return;
150 G4double E2 = kinEnergyProd * 1.000001; << 98 }
151 G4double dE = (E2 - E1); << 99
152 G4double sigma1, sigma2; << 100 //Sample secondary energy
153 fDirectModel = fBraggIonDirectEMModel; << 101 //-----------------------
154 if(kinEnergyProjScaled > 2. * MeV && !fUse << 102 G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, IsScatProjToProjCase);
155 fDirectModel = fBetheBlochDirectEMModel; << 103 CorrectPostStepWeight(fParticleChange,
156 sigma1 = fDirectModel->ComputeCrossSection << 104 aTrack.GetWeight(),
157 fDirectPrimaryPart, kinEnergyProj, Z, A, << 105 adjointPrimKinEnergy,
158 sigma2 = fDirectModel->ComputeCrossSection << 106 projectileKinEnergy,
159 fDirectPrimaryPart, kinEnergyProj, Z, A, << 107 IsScatProjToProjCase); //Caution !!!this weight correction should be always applied
160 << 108
161 dSigmadEprod = (sigma1 - sigma2) / dE; << 109
162 << 110 //Kinematic:
163 if(dSigmadEprod > 1.) << 111 //we consider a two body elastic scattering for the forward processes where the projectile knock on an e- at rest and gives
164 { << 112 // him part of its energy
165 G4cout << "sigma1 " << kinEnergyProj / M << 113 //----------------------------------------------------------------------------------------
166 << '\t' << sigma1 << G4endl; << 114
167 G4cout << "sigma2 " << kinEnergyProj / M << 115 G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
168 << '\t' << sigma2 << G4endl; << 116 G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
169 G4cout << "dsigma " << kinEnergyProj / M << 117 G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
170 << '\t' << dSigmadEprod << G4endl << 118
171 } << 119
>> 120
>> 121 //Companion
>> 122 //-----------
>> 123 G4double companionM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
>> 124 if (IsScatProjToProjCase) {
>> 125 companionM0=theAdjEquivOfDirectSecondPartDef->GetPDGMass();
>> 126 }
>> 127 G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
>> 128 G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;
>> 129
>> 130
>> 131 //Projectile momentum
>> 132 //--------------------
>> 133 G4double P_parallel = (adjointPrimP*adjointPrimP + projectileP2 - companionP2)/(2.*adjointPrimP);
>> 134 G4double P_perp = std::sqrt( projectileP2 - P_parallel*P_parallel);
>> 135 G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
>> 136 G4double phi =G4UniformRand()*2.*3.1415926;
>> 137 G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
>> 138 projectileMomentum.rotateUz(dir_parallel);
>> 139
>> 140
>> 141
>> 142 if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
>> 143 fParticleChange->ProposeTrackStatus(fStopAndKill);
>> 144 fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
>> 145 //G4cout<<"projectileMomentum "<<projectileMomentum<<G4endl;
>> 146 }
>> 147 else {
>> 148 fParticleChange->ProposeEnergy(projectileKinEnergy);
>> 149 fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
>> 150 }
>> 151
172 152
173 if(fDirectModel == fBetheBlochDirectEMMode << 153
174 { <<
175 // correction of differential cross sect <<
176 // the suppression of particle at second <<
177 // Bethe Bloch Model. This correction co <<
178 // probability function used to test the <<
179 // code taken from G4BetheBlochModel::Sa <<
180 <<
181 G4double deltaKinEnergy = kinEnergyProd; <<
182 <<
183 G4double x = fFormFact * deltaKinEnergy; <<
184 if(x > 1.e-6) <<
185 { <<
186 G4double totEnergy = kinEnergyProj + f <<
187 G4double etot2 = totEnergy * totEn <<
188 G4double beta2 = kinEnergyProj * (kinE <<
189 G4double f1 = 0.0; <<
190 G4double f = 1.0 - beta2 * deltaKi <<
191 if(0.5 == fSpin) <<
192 { <<
193 f1 = 0.5 * deltaKinEnergy * deltaKin <<
194 f += f1; <<
195 } <<
196 G4double x1 = 1.0 + x; <<
197 G4double gg = 1.0 / (x1 * x1); <<
198 if(0.5 == fSpin) <<
199 { <<
200 G4double x2 = <<
201 0.5 * electron_mass_c2 * deltaKinE <<
202 gg *= (1.0 + fMagMoment2 * (x2 - f1 <<
203 } <<
204 if(gg > 1.0) <<
205 { <<
206 G4cout << "### G4BetheBlochModel in <<
207 << G4endl; <<
208 gg = 1.; <<
209 } <<
210 dSigmadEprod *= gg; <<
211 } <<
212 } <<
213 } <<
214 return dSigmadEprod; <<
215 } <<
216 154
>> 155 }
>> 156
217 ////////////////////////////////////////////// 157 ////////////////////////////////////////////////////////////////////////////////
218 void G4AdjointIonIonisationModel::SetIon(G4Par << 158 //
219 G4Par << 159 G4double G4AdjointIonIonisationModel::DiffCrossSectionPerAtomPrimToSecond(
220 { << 160 G4double kinEnergyProj,
221 fDirectPrimaryPart = fwd_ion; << 161 G4double kinEnergyProd,
222 fAdjEquivDirectPrimPart = adj_ion; << 162 G4double Z,
>> 163 G4double A)
>> 164 {//Probably that here the Bragg Model should be also used for kinEnergyProj/nuc<2MeV
>> 165
>> 166
>> 167
>> 168 G4double dSigmadEprod=0;
>> 169 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
>> 170 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
>> 171
>> 172 G4double kinEnergyProjScaled = massRatio*kinEnergyProj;
>> 173
>> 174
>> 175 if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ //the produced particle should have a kinetic energy smaller than the projectile
>> 176 G4double Tmax=kinEnergyProj;
>> 177
>> 178 G4double E1=kinEnergyProd;
>> 179 G4double E2=kinEnergyProd*1.000001;
>> 180 G4double dE=(E2-E1);
>> 181 G4double sigma1,sigma2;
>> 182 theDirectEMModel =theBraggIonDirectEMModel;
>> 183 if (kinEnergyProjScaled >2.*MeV && !use_only_bragg) theDirectEMModel = theBetheBlochDirectEMModel; //Bethe Bloch Model
>> 184 sigma1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
>> 185 sigma2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
>> 186
>> 187 dSigmadEprod=(sigma1-sigma2)/dE;
>> 188
>> 189 //G4double chargeSqRatio = currentModel->GetChargeSquareRatio(theDirectPrimaryPartDef,currentMaterial,E);
>> 190
>> 191
>> 192
>> 193 if (dSigmadEprod>1.) {
>> 194 G4cout<<"sigma1 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma1<<G4endl;
>> 195 G4cout<<"sigma2 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma2<<G4endl;
>> 196 G4cout<<"dsigma "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<dSigmadEprod<<G4endl;
>> 197
>> 198 }
>> 199
>> 200
>> 201
>> 202
>> 203
>> 204
>> 205 if (theDirectEMModel == theBetheBlochDirectEMModel ){
>> 206 //correction of differential cross section at high energy to correct for the suppression of particle at secondary at high
>> 207 //energy used in the Bethe Bloch Model. This correction consist to multiply by g the probability function used
>> 208 //to test the rejection of a secondary
>> 209 //-------------------------
>> 210
>> 211 //Source code taken from G4BetheBlochModel::SampleSecondaries
>> 212
>> 213 G4double deltaKinEnergy = kinEnergyProd;
>> 214
>> 215 //Part of the taken code
>> 216 //----------------------
>> 217
>> 218
>> 219
>> 220 // projectile formfactor - suppresion of high energy
>> 221 // delta-electron production at high energy
>> 222
>> 223
>> 224 G4double x = formfact*deltaKinEnergy;
>> 225 if(x > 1.e-6) {
>> 226 G4double totEnergy = kinEnergyProj + mass;
>> 227 G4double etot2 = totEnergy*totEnergy;
>> 228 G4double beta2 = kinEnergyProj*(kinEnergyProj + 2.0*mass)/etot2;
>> 229 G4double f;
>> 230 G4double f1 = 0.0;
>> 231 f = 1.0 - beta2*deltaKinEnergy/Tmax;
>> 232 if( 0.5 == spin ) {
>> 233 f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
>> 234 f += f1;
>> 235 }
>> 236 G4double x1 = 1.0 + x;
>> 237 G4double g = 1.0/(x1*x1);
>> 238 if( 0.5 == spin ) {
>> 239 G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
>> 240 g *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
>> 241 }
>> 242 if(g > 1.0) {
>> 243 G4cout << "### G4BetheBlochModel in Adjoint Sim WARNING: g= " << g
>> 244 << G4endl;
>> 245 g=1.;
>> 246 }
>> 247 //G4cout<<"g"<<g<<G4endl;
>> 248 dSigmadEprod*=g;
>> 249 }
>> 250 }
>> 251
>> 252 }
223 253
>> 254 return dSigmadEprod;
>> 255 }
>> 256 //////////////////////////////////////////////////////////////////////////////////////////////
>> 257 //
>> 258 void G4AdjointIonIonisationModel::SetIon(G4ParticleDefinition* adj_ion, G4ParticleDefinition* fwd_ion)
>> 259 { theDirectPrimaryPartDef =fwd_ion;
>> 260 theAdjEquivOfDirectPrimPartDef =adj_ion;
>> 261
224 DefineProjectileProperty(); 262 DefineProjectileProperty();
225 } 263 }
226 << 264 //////////////////////////////////////////////////////////////////////////////////////////////
227 ////////////////////////////////////////////// << 265 //
228 void G4AdjointIonIonisationModel::CorrectPostS << 266 void G4AdjointIonIonisationModel::CorrectPostStepWeight(G4ParticleChange* fParticleChange, G4double old_weight,
229 G4ParticleChange* fParticleChange, G4double << 267 G4double adjointPrimKinEnergy, G4double projectileKinEnergy,G4bool )
230 G4double adjointPrimKinEnergy, G4double proj <<
231 { 268 {
232 // It is needed because the direct cross sec << 269 //It is needed because the direct cross section used to compute the differential cross section is not the one used in
233 // differential cross section is not the one << 270 // the direct model where the GenericIon stuff is considered with correction of effective charge. In the direct model the samnepl of secondaries does
234 // the direct model where the GenericIon stu << 271 // not reflect the integral cross section. The integral fwd cross section that we used to compute the differential CS
235 // of effective charge. In the direct model << 272 // match the sample of secondaries in the forward case despite the fact that its is not the same total CS than in the FWD case. For this reasion an extra
236 // not reflect the integral cross section. T << 273 // weight correction is needed at the end.
237 // we used to compute the differential CS ma << 274
238 // the forward case despite the fact that it << 275
239 // the FWD case. For this reason an extra we << 276 G4double new_weight=old_weight;
240 // end. << 277
241 << 278 //the correction of CS due to the problem explained above
242 G4double new_weight = old_weight; << 279 G4double kinEnergyProjScaled = massRatio*projectileKinEnergy;
243 << 280 theDirectEMModel =theBraggIonDirectEMModel;
244 // the correction of CS due to the problem e << 281 if (kinEnergyProjScaled >2.*MeV && !use_only_bragg) theDirectEMModel = theBetheBlochDirectEMModel; //Bethe Bloch Model
245 G4double kinEnergyProjScaled = fMassRatio * << 282 G4double UsedFwdCS=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,projectileKinEnergy,1,1 ,currentTcutForDirectSecond,1.e20);
246 fDirectModel = fBraggIonDire << 283 G4double chargeSqRatio =1.;
247 if(kinEnergyProjScaled > 2. * MeV && !fUseOn << 284 if (chargeSquare>1.) chargeSqRatio = theDirectEMModel->GetChargeSquareRatio(theDirectPrimaryPartDef,currentMaterial,projectileKinEnergy);
248 fDirectModel = fBetheBlochDirectEMModel; << 285 G4double CorrectFwdCS = chargeSqRatio*theDirectEMModel->ComputeCrossSectionPerAtom(G4GenericIon::GenericIon(),kinEnergyProjScaled,1,1 ,currentTcutForDirectSecond,1.e20);
249 G4double UsedFwdCS = fDirectModel->ComputeCr << 286 if (UsedFwdCS >0) new_weight*= CorrectFwdCS/UsedFwdCS;//May be some check is needed if UsedFwdCS ==0 probably that then we should avoid a secondary to be produced,
250 fDirectPrimaryPart, projectileKinEnergy, 1 << 287
251 G4double chargeSqRatio = 1.; << 288
252 if(fChargeSquare > 1.) << 289 //additional CS crorrection needed for cross section biasing in general.
253 chargeSqRatio = fDirectModel->GetChargeSqu << 290 //May be wrong for ions!!! Most of the time not used!
254 fDirectPrimaryPart, fCurrentMaterial, pr << 291 G4double w_corr =1./CS_biasing_factor;
255 G4double CorrectFwdCS = << 292 w_corr*=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
256 chargeSqRatio * fDirectModel->ComputeCross << 293 new_weight*=w_corr;
257 G4GenericIon::GenericIon << 294
258 fTcutSecond, 1.e20); << 295 new_weight*=projectileKinEnergy/adjointPrimKinEnergy;
259 // May be some check is needed if UsedFwdCS << 296
260 // avoid a secondary to be produced, << 297 fParticleChange->SetParentWeightByProcess(false);
261 if(UsedFwdCS > 0.) << 298 fParticleChange->SetSecondaryWeightByProcess(false);
262 new_weight *= CorrectFwdCS / UsedFwdCS; << 299 fParticleChange->ProposeParentWeight(new_weight);
263 <<
264 // additional CS correction needed for cross <<
265 // May be wrong for ions. Most of the time n <<
266 new_weight *= <<
267 G4AdjointCSManager::GetAdjointCSManager()- <<
268 fCsBiasingFactor; <<
269 <<
270 new_weight *= projectileKinEnergy / adjointP <<
271 <<
272 fParticleChange->SetParentWeightByProcess(fa <<
273 fParticleChange->SetSecondaryWeightByProcess <<
274 fParticleChange->ProposeParentWeight(new_wei <<
275 } 300 }
276 301
277 ////////////////////////////////////////////// <<
278 void G4AdjointIonIonisationModel::DefineProjec <<
279 { <<
280 // Slightly modified code taken from G4Bethe <<
281 G4String pname = fDirectPrimaryPart->GetPart <<
282 302
283 fMass = fDirectPrimaryPart->GetPDG << 303 //////////////////////////////////////////////////////////////////////////////////////////////
284 fMassRatio = G4GenericIon::GenericIon() << 304 //
285 fSpin = fDirectPrimaryPart->GetPDG << 305 void G4AdjointIonIonisationModel::DefineProjectileProperty()
286 G4double q = fDirectPrimaryPart->GetPDG << 306 {
287 fChargeSquare = q * q; << 307 //Slightly modified code taken from G4BetheBlochModel::SetParticle
288 fRatio = electron_mass_c2 / fMass; << 308 //------------------------------------------------
289 fOnePlusRatio2 = (1. + fRatio) * (1. + fRat << 309 G4String pname = theDirectPrimaryPartDef->GetParticleName();
290 fOneMinusRatio2 = (1. - fRatio) * (1. - fRat << 310 if (theDirectPrimaryPartDef->GetParticleType() == "nucleus" &&
291 G4double magmom = fDirectPrimaryPart->GetPDG << 311 pname != "deuteron" && pname != "triton") {
292 (0.5 * eplus * hbar_Planck << 312 isIon = true;
293 fMagMoment2 = magmom * magmom - 1.0; <<
294 if(fDirectPrimaryPart->GetLeptonNumber() == <<
295 { <<
296 G4double x = 0.8426 * GeV; <<
297 if(fSpin == 0.0 && fMass < GeV) <<
298 { <<
299 x = 0.736 * GeV; <<
300 } <<
301 else if(fMass > GeV) <<
302 { <<
303 x /= G4NistManager::Instance()->GetZ13(f <<
304 } 313 }
305 fFormFact = 2.0 * electron_mass_c2 / (x * << 314
306 } << 315 mass = theDirectPrimaryPartDef->GetPDGMass();
>> 316 massRatio= G4GenericIon::GenericIon()->GetPDGMass()/mass;
>> 317 mass_ratio_projectile = massRatio;
>> 318 spin = theDirectPrimaryPartDef->GetPDGSpin();
>> 319 G4double q = theDirectPrimaryPartDef->GetPDGCharge()/eplus;
>> 320 chargeSquare = q*q;
>> 321 ratio = electron_mass_c2/mass;
>> 322 ratio2 = ratio*ratio;
>> 323 one_plus_ratio_2=(1+ratio)*(1+ratio);
>> 324 one_minus_ratio_2=(1-ratio)*(1-ratio);
>> 325 G4double magmom = theDirectPrimaryPartDef->GetPDGMagneticMoment()
>> 326 *mass/(0.5*eplus*hbar_Planck*c_squared);
>> 327 magMoment2 = magmom*magmom - 1.0;
>> 328 formfact = 0.0;
>> 329 if(theDirectPrimaryPartDef->GetLeptonNumber() == 0) {
>> 330 G4double x = 0.8426*GeV;
>> 331 if(spin == 0.0 && mass < GeV) {x = 0.736*GeV;}
>> 332 else if(mass > GeV) {
>> 333 x /= G4NistManager::Instance()->GetZ13(mass/proton_mass_c2);
>> 334 // tlimit = 51.2*GeV*A13[iz]*A13[iz];
>> 335 }
>> 336 formfact = 2.0*electron_mass_c2/(x*x);
>> 337 tlimit = 2.0/formfact;
>> 338 }
307 } 339 }
308 340
>> 341
309 ////////////////////////////////////////////// 342 //////////////////////////////////////////////////////////////////////////////
310 G4double G4AdjointIonIonisationModel::GetSecon << 343 //
311 G4double primAdjEnergy) << 344 G4double G4AdjointIonIonisationModel::GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
312 { << 345 {
313 return primAdjEnergy * fOnePlusRatio2 / << 346 G4double Tmax=PrimAdjEnergy*one_plus_ratio_2/(one_minus_ratio_2-2.*ratio*PrimAdjEnergy/mass);
314 (fOneMinusRatio2 - 2. * fRatio * prim << 347 return Tmax;
315 } 348 }
316 <<
317 ////////////////////////////////////////////// 349 //////////////////////////////////////////////////////////////////////////////
318 G4double G4AdjointIonIonisationModel::GetSecon << 350 //
319 G4double primAdjEnergy, G4double tcut) << 351 G4double G4AdjointIonIonisationModel::GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy,G4double Tcut)
320 { << 352 { return PrimAdjEnergy+Tcut;
321 return primAdjEnergy + tcut; <<
322 } 353 }
323 <<
324 ////////////////////////////////////////////// 354 //////////////////////////////////////////////////////////////////////////////
325 G4double G4AdjointIonIonisationModel::GetSecon << 355 //
326 G4double) << 356 G4double G4AdjointIonIonisationModel::GetSecondAdjEnergyMaxForProdToProjCase(G4double )
327 { << 357 { return HighEnergyLimit;
328 return GetHighEnergyLimit(); <<
329 } 358 }
330 <<
331 ////////////////////////////////////////////// 359 //////////////////////////////////////////////////////////////////////////////
332 G4double G4AdjointIonIonisationModel::GetSecon << 360 //
333 G4double primAdjEnergy) << 361 G4double G4AdjointIonIonisationModel::GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
334 { << 362 { G4double Tmin= (2*PrimAdjEnergy-4*mass + std::sqrt(4.*PrimAdjEnergy*PrimAdjEnergy +16.*mass*mass + 8.*PrimAdjEnergy*mass*(1/ratio +ratio)))/4.;
335 return (2. * primAdjEnergy - 4. * fMass + << 363 return Tmin;
336 std::sqrt(4. * primAdjEnergy * primA <<
337 8. * primAdjEnergy * fMass <<
338 4.; <<
339 } 364 }
340 365