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