| 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 "G4AdjointhIonisationModel.hh" 27 #include "G4AdjointhIonisationModel.hh"
28 28
>> 29 #include "G4PhysicalConstants.hh"
>> 30 #include "G4SystemOfUnits.hh"
29 #include "G4AdjointCSManager.hh" 31 #include "G4AdjointCSManager.hh"
>> 32 #include "G4Integrator.hh"
>> 33 #include "G4TrackStatus.hh"
>> 34 #include "G4ParticleChange.hh"
30 #include "G4AdjointElectron.hh" 35 #include "G4AdjointElectron.hh"
31 #include "G4AdjointProton.hh" 36 #include "G4AdjointProton.hh"
>> 37 #include "G4AdjointInterpolator.hh"
32 #include "G4BetheBlochModel.hh" 38 #include "G4BetheBlochModel.hh"
33 #include "G4BraggModel.hh" 39 #include "G4BraggModel.hh"
34 #include "G4NistManager.hh" <<
35 #include "G4ParticleChange.hh" <<
36 #include "G4PhysicalConstants.hh" <<
37 #include "G4Proton.hh" 40 #include "G4Proton.hh"
38 #include "G4SystemOfUnits.hh" << 41 #include "G4NistManager.hh"
39 #include "G4TrackStatus.hh" <<
40 42
41 ////////////////////////////////////////////// 43 ////////////////////////////////////////////////////////////////////////////////
42 G4AdjointhIonisationModel::G4AdjointhIonisatio << 44 //
43 : G4VEmAdjointModel("Adjoint_hIonisation") << 45 G4AdjointhIonisationModel::G4AdjointhIonisationModel(G4ParticleDefinition* projectileDefinition):
44 { << 46 G4VEmAdjointModel("Adjoint_hIonisation")
45 fUseMatrix = true; << 47 {
46 fUseMatrixPerElement = true; << 48
47 fApplyCutInRange = true; << 49
48 fOneMatrixForAllElements = true; <<
49 fSecondPartSameType = false; <<
50 <<
51 // The direct EM Model is taken as BetheBloc <<
52 // computation of the differential cross sec <<
53 // The Bragg model could be used as an alter <<
54 // differential cross section <<
55 <<
56 fDirectModel = new G4BetheBloch <<
57 fBraggDirectEMModel = new G4BraggModel <<
58 fAdjEquivDirectSecondPart = G4AdjointElectro <<
59 fDirectPrimaryPart = pDef; <<
60 <<
61 if(pDef == G4Proton::Proton()) <<
62 { <<
63 fAdjEquivDirectPrimPart = G4AdjointProton: <<
64 } <<
65 50
>> 51 UseMatrix =true;
>> 52 UseMatrixPerElement = true;
>> 53 ApplyCutInRange = true;
>> 54 UseOnlyOneMatrixForAllElements = true;
>> 55 CS_biasing_factor =1.;
>> 56 second_part_of_same_type =false;
>> 57
>> 58 //The direct EM Modfel is taken has BetheBloch it is only used for the computation
>> 59 // of the differential cross section.
>> 60 //The Bragg model could be used as an alternative as it offers the same differential cross section
>> 61
>> 62 theDirectEMModel = new G4BetheBlochModel(projectileDefinition);
>> 63 theBraggDirectEMModel = new G4BraggModel(projectileDefinition);
>> 64 theAdjEquivOfDirectSecondPartDef=G4AdjointElectron::AdjointElectron();
>> 65
>> 66 theDirectPrimaryPartDef = projectileDefinition;
>> 67 theAdjEquivOfDirectPrimPartDef = 0;
>> 68 if (projectileDefinition == G4Proton::Proton()) {
>> 69 theAdjEquivOfDirectPrimPartDef = G4AdjointProton::AdjointProton();
>> 70
>> 71 }
>> 72
66 DefineProjectileProperty(); 73 DefineProjectileProperty();
67 } 74 }
68 <<
69 ////////////////////////////////////////////// 75 ////////////////////////////////////////////////////////////////////////////////
70 G4AdjointhIonisationModel::~G4AdjointhIonisati << 76 //
>> 77 G4AdjointhIonisationModel::~G4AdjointhIonisationModel()
>> 78 {;}
>> 79
71 80
72 ////////////////////////////////////////////// 81 ////////////////////////////////////////////////////////////////////////////////
73 void G4AdjointhIonisationModel::SampleSecondar << 82 //
74 const G4Track& aTrack, G4bool isScatProjToPr << 83 void G4AdjointhIonisationModel::SampleSecondaries(const G4Track& aTrack,
75 G4ParticleChange* fParticleChange) << 84 G4bool IsScatProjToProjCase,
76 { << 85 G4ParticleChange* fParticleChange)
77 if(!fUseMatrix) << 86 {
78 return RapidSampleSecondaries(aTrack, isSc << 87 if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange);
79 << 88
80 const G4DynamicParticle* theAdjointPrimary = << 89 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
81 << 90
82 // Elastic inverse scattering << 91 //Elastic inverse scattering
83 G4double adjointPrimKinEnergy = theAdjointPr << 92 //---------------------------------------------------------
84 G4double adjointPrimP = theAdjointPr << 93 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
85 << 94 G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
86 if(adjointPrimKinEnergy > GetHighEnergyLimit << 95
87 { << 96 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
88 return; << 97 return;
89 } << 98 }
>> 99
>> 100 //Sample secondary energy
>> 101 //-----------------------
>> 102 G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, IsScatProjToProjCase);
>> 103 CorrectPostStepWeight(fParticleChange,
>> 104 aTrack.GetWeight(),
>> 105 adjointPrimKinEnergy,
>> 106 projectileKinEnergy,
>> 107 IsScatProjToProjCase); //Caution !!!this weight correction should be always applied
>> 108
>> 109
>> 110 //Kinematic:
>> 111 //we consider a two body elastic scattering for the forward processes where the projectile knock on an e- at rest and gives
>> 112 // him part of its energy
>> 113 //----------------------------------------------------------------------------------------
>> 114
>> 115 G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
>> 116 G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
>> 117 G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
>> 118
>> 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
>> 152
>> 153
90 154
91 // Sample secondary energy <<
92 G4double projectileKinEnergy = <<
93 SampleAdjSecEnergyFromCSMatrix(adjointPrim <<
94 CorrectPostStepWeight(fParticleChange, aTrac <<
95 adjointPrimKinEnergy, <<
96 isScatProjToProj); <<
97 // Caution!!! this weight correction should <<
98 <<
99 // Kinematic: <<
100 // we consider a two body elastic scattering <<
101 // the projectile knock on an e- at rest and <<
102 G4double projectileM0 = fAdjEquivDi <<
103 G4double projectileTotalEnergy = projectileM <<
104 G4double projectileP2 = <<
105 projectileTotalEnergy * projectileTotalEne <<
106 <<
107 // Companion <<
108 G4double companionM0 = fAdjEquivDirectPrimPa <<
109 if(isScatProjToProj) <<
110 { <<
111 companionM0 = fAdjEquivDirectSecondPart->G <<
112 } <<
113 G4double companionTotalEnergy = <<
114 companionM0 + projectileKinEnergy - adjoin <<
115 G4double companionP2 = <<
116 companionTotalEnergy * companionTotalEnerg <<
117 <<
118 // Projectile momentum <<
119 G4double P_parallel = <<
120 (adjointPrimP * adjointPrimP + projectileP <<
121 (2. * adjointPrimP); <<
122 G4double P_perp = std::sqrt(projectileP2 - P <<
123 G4ThreeVector dir_parallel = theAdjointPrima <<
124 G4double phi = G4UniformRand() <<
125 G4ThreeVector projectileMomentum = <<
126 G4ThreeVector(P_perp * std::cos(phi), P_pe <<
127 projectileMomentum.rotateUz(dir_parallel); <<
128 <<
129 if(!isScatProjToProj) <<
130 { // kill the primary and add a secondary <<
131 fParticleChange->ProposeTrackStatus(fStopA <<
132 fParticleChange->AddSecondary( <<
133 new G4DynamicParticle(fAdjEquivDirectPri <<
134 } <<
135 else <<
136 { <<
137 fParticleChange->ProposeEnergy(projectileK <<
138 fParticleChange->ProposeMomentumDirection( <<
139 } <<
140 } 155 }
141 156
142 ////////////////////////////////////////////// 157 ////////////////////////////////////////////////////////////////////////////////
143 void G4AdjointhIonisationModel::RapidSampleSec << 158 //
144 const G4Track& aTrack, G4bool isScatProjToPr << 159 void G4AdjointhIonisationModel::RapidSampleSecondaries(const G4Track& aTrack,
145 G4ParticleChange* fParticleChange) << 160 G4bool IsScatProjToProjCase,
146 { << 161 G4ParticleChange* fParticleChange)
147 const G4DynamicParticle* theAdjointPrimary = << 162 {
148 DefineCurrentMaterial(aTrack.GetMaterialCuts << 163
149 << 164 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
150 G4double adjointPrimKinEnergy = theAdjointPr << 165 DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
151 G4double adjointPrimP = theAdjointPr << 166
152 << 167
153 if(adjointPrimKinEnergy > GetHighEnergyLimit << 168 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
154 { << 169 G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
155 return; << 170
156 } << 171 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
>> 172 return;
>> 173 }
>> 174
>> 175 G4double projectileKinEnergy =0.;
>> 176 G4double eEnergy=0.;
>> 177 G4double newCS=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2*mass;
>> 178 if (!IsScatProjToProjCase){//1/E^2 distribution
>> 179
>> 180 eEnergy=adjointPrimKinEnergy;
>> 181 G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
>> 182 G4double Emin= GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);
>> 183 if (Emin>=Emax) return;
>> 184 G4double a=1./Emax;
>> 185 G4double b=1./Emin;
>> 186 newCS=newCS*(b-a)/eEnergy;
>> 187
>> 188 projectileKinEnergy =1./(b- (b-a)*G4UniformRand());
>> 189
>> 190
>> 191 }
>> 192 else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
>> 193 G4double Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond);
>> 194 if (Emin>=Emax) return;
>> 195 G4double diff1=Emin-adjointPrimKinEnergy;
>> 196 G4double diff2=Emax-adjointPrimKinEnergy;
>> 197
>> 198 G4double t1=adjointPrimKinEnergy*(1./diff1-1./diff2);
>> 199 G4double t2=adjointPrimKinEnergy*(1./Emin-1./Emax);
>> 200 /*G4double f31=diff1/Emin;
>> 201 G4double f32=diff2/Emax/f31;*/
>> 202 G4double t3=2.*std::log(Emax/Emin);
>> 203 G4double sum_t=t1+t2+t3;
>> 204 newCS=newCS*sum_t/adjointPrimKinEnergy/adjointPrimKinEnergy;
>> 205 G4double t=G4UniformRand()*sum_t;
>> 206 if (t <=t1 ){
>> 207 G4double q= G4UniformRand()*t1/adjointPrimKinEnergy ;
>> 208 projectileKinEnergy =adjointPrimKinEnergy +1./(1./diff1-q);
>> 209
>> 210 }
>> 211 else if (t <=t2 ) {
>> 212 G4double q= G4UniformRand()*t2/adjointPrimKinEnergy;
>> 213 projectileKinEnergy =1./(1./Emin-q);
>> 214 }
>> 215 else {
>> 216 projectileKinEnergy=Emin*std::pow(Emax/Emin,G4UniformRand());
>> 217
>> 218 }
>> 219 eEnergy=projectileKinEnergy-adjointPrimKinEnergy;
>> 220
>> 221
>> 222 }
>> 223
>> 224
>> 225
>> 226 G4double diffCS_perAtom_Used=twopi_mc2_rcl2*mass*adjointPrimKinEnergy/projectileKinEnergy/projectileKinEnergy/eEnergy/eEnergy;
>> 227
>> 228
>> 229
>> 230 //Weight correction
>> 231 //-----------------------
>> 232 //First w_corr is set to the ratio between adjoint total CS and fwd total CS
>> 233 G4double w_corr=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
>> 234
>> 235 //G4cout<<w_corr<<G4endl;
>> 236 w_corr*=newCS/lastCS;
>> 237 //G4cout<<w_corr<<G4endl;
>> 238 //Then another correction is needed due to the fact that a biaised differential CS has been used rather than the one consistent with the direct model
>> 239 //Here we consider the true diffCS as the one obtained by the numerical differentiation over Tcut of the direct CS
>> 240
>> 241 G4double diffCS = DiffCrossSectionPerAtomPrimToSecond(projectileKinEnergy, eEnergy,1,1);
>> 242 w_corr*=diffCS/diffCS_perAtom_Used;
>> 243 //G4cout<<w_corr<<G4endl;
>> 244
>> 245 G4double new_weight = aTrack.GetWeight()*w_corr;
>> 246 fParticleChange->SetParentWeightByProcess(false);
>> 247 fParticleChange->SetSecondaryWeightByProcess(false);
>> 248 fParticleChange->ProposeParentWeight(new_weight);
>> 249
>> 250
>> 251
>> 252
>> 253 //Kinematic:
>> 254 //we consider a two body elastic scattering for the forward processes where the projectile knock on an e- at rest and gives
>> 255 // him part of its energy
>> 256 //----------------------------------------------------------------------------------------
>> 257
>> 258 G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
>> 259 G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
>> 260 G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
>> 261
>> 262
>> 263
>> 264 //Companion
>> 265 //-----------
>> 266 G4double companionM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
>> 267 if (IsScatProjToProjCase) {
>> 268 companionM0=theAdjEquivOfDirectSecondPartDef->GetPDGMass();
>> 269 }
>> 270 G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
>> 271 G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;
>> 272
>> 273
>> 274 //Projectile momentum
>> 275 //--------------------
>> 276 G4double P_parallel = (adjointPrimP*adjointPrimP + projectileP2 - companionP2)/(2.*adjointPrimP);
>> 277 G4double P_perp = std::sqrt( projectileP2 - P_parallel*P_parallel);
>> 278 G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
>> 279 G4double phi =G4UniformRand()*2.*3.1415926;
>> 280 G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
>> 281 projectileMomentum.rotateUz(dir_parallel);
>> 282
>> 283
>> 284
>> 285 if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
>> 286 fParticleChange->ProposeTrackStatus(fStopAndKill);
>> 287 fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
>> 288 //G4cout<<"projectileMomentum "<<projectileMomentum<<G4endl;
>> 289 }
>> 290 else {
>> 291 fParticleChange->ProposeEnergy(projectileKinEnergy);
>> 292 fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
>> 293 }
>> 294
157 295
158 G4double projectileKinEnergy = 0.; << 296
159 G4double eEnergy = 0.; << 297
160 G4double newCS = <<
161 fCurrentMaterial->GetElectronDensity() * t <<
162 if(!isScatProjToProj) <<
163 { // 1/E^2 distribution <<
164 <<
165 eEnergy = adjointPrimKinEnergy; <<
166 G4double Emax = GetSecondAdjEnergyMaxForPr <<
167 G4double Emin = GetSecondAdjEnergyMinForPr <<
168 if(Emin >= Emax) <<
169 return; <<
170 G4double a = 1. / Emax; <<
171 G4double b = 1. / Emin; <<
172 newCS = newCS * (b - a) / eEnergy; <<
173 298
174 projectileKinEnergy = 1. / (b - (b - a) * <<
175 } <<
176 else <<
177 { <<
178 G4double Emax = <<
179 GetSecondAdjEnergyMaxForScatProjToProj(a <<
180 G4double Emin = <<
181 GetSecondAdjEnergyMinForScatProjToProj(a <<
182 if(Emin >= Emax) <<
183 return; <<
184 G4double diff1 = Emin - adjointPrimKinEner <<
185 G4double diff2 = Emax - adjointPrimKinEner <<
186 <<
187 G4double t1 = adjointPrimKinEnergy * (1 <<
188 G4double t2 = adjointPrimKinEnergy * (1 <<
189 G4double t3 = 2. * std::log(Emax / Emin <<
190 G4double sum_t = t1 + t2 + t3; <<
191 newCS = newCS * sum_t / adjointPrimKi <<
192 G4double t = G4UniformRand() * sum_t; <<
193 if(t <= t1) <<
194 { <<
195 G4double q = G4UniformRand() * <<
196 projectileKinEnergy = adjointPrimKinEner <<
197 } <<
198 else if(t <= t2) <<
199 { <<
200 G4double q = G4UniformRand() * <<
201 projectileKinEnergy = 1. / (1. / Emin - <<
202 } <<
203 else <<
204 { <<
205 projectileKinEnergy = Emin * std::pow(Em <<
206 } <<
207 eEnergy = projectileKinEnergy - adjointPri <<
208 } <<
209 299
210 G4double diffCS_perAtom_Used = twopi_mc2_rcl <<
211 projectileKin <<
212 eEnergy / eEn <<
213 <<
214 // Weight correction <<
215 // First w_corr is set to the ratio between <<
216 G4double w_corr = <<
217 G4AdjointCSManager::GetAdjointCSManager()- <<
218 <<
219 w_corr *= newCS / fLastCS; <<
220 // Then another correction is needed due to <<
221 // differential CS has been used rather than <<
222 // direct model. Here we consider the true d <<
223 // numerical differentiation over Tcut of th <<
224 <<
225 G4double diffCS = <<
226 DiffCrossSectionPerAtomPrimToSecond(projec <<
227 w_corr *= diffCS / diffCS_perAtom_Used; <<
228 <<
229 if (isScatProjToProj && fTcutSecond>0.005) w <<
230 <<
231 G4double new_weight = aTrack.GetWeight() * w <<
232 fParticleChange->SetParentWeightByProcess(fa <<
233 fParticleChange->SetSecondaryWeightByProcess <<
234 fParticleChange->ProposeParentWeight(new_wei <<
235 <<
236 // Kinematic: <<
237 // we consider a two body elastic scattering <<
238 // the projectile knocks on an e- at rest an <<
239 G4double projectileM0 = fAdjEquivDi <<
240 G4double projectileTotalEnergy = projectileM <<
241 G4double projectileP2 = <<
242 projectileTotalEnergy * projectileTotalEne <<
243 <<
244 // Companion <<
245 G4double companionM0 = fAdjEquivDirectPrimPa <<
246 if(isScatProjToProj) <<
247 { <<
248 companionM0 = fAdjEquivDirectSecondPart->G <<
249 } <<
250 G4double companionTotalEnergy = <<
251 companionM0 + projectileKinEnergy - adjoin <<
252 G4double companionP2 = <<
253 companionTotalEnergy * companionTotalEnerg <<
254 <<
255 // Projectile momentum <<
256 G4double P_parallel = <<
257 (adjointPrimP * adjointPrimP + projectileP <<
258 (2. * adjointPrimP); <<
259 G4double P_perp = std::sqrt(projectileP2 - P <<
260 G4ThreeVector dir_parallel = theAdjointPrima <<
261 G4double phi = G4UniformRand() <<
262 G4ThreeVector projectileMomentum = <<
263 G4ThreeVector(P_perp * std::cos(phi), P_pe <<
264 projectileMomentum.rotateUz(dir_parallel); <<
265 <<
266 if(!isScatProjToProj) <<
267 { // kill the primary and add a secondary <<
268 fParticleChange->ProposeTrackStatus(fStopA <<
269 fParticleChange->AddSecondary( <<
270 new G4DynamicParticle(fAdjEquivDirectPri <<
271 } <<
272 else <<
273 { <<
274 fParticleChange->ProposeEnergy(projectileK <<
275 fParticleChange->ProposeMomentumDirection( <<
276 } <<
277 } <<
278 300
>> 301 }
>> 302
279 ////////////////////////////////////////////// 303 ////////////////////////////////////////////////////////////////////////////////
>> 304 //
280 G4double G4AdjointhIonisationModel::DiffCrossS 305 G4double G4AdjointhIonisationModel::DiffCrossSectionPerAtomPrimToSecond(
281 G4double kinEnergyProj, G4double kinEnergyPr << 306 G4double kinEnergyProj,
282 { // Probably here the Bragg Model should be << 307 G4double kinEnergyProd,
283 // kinEnergyProj/nuc < 2 MeV << 308 G4double Z,
284 << 309 G4double A)
285 G4double dSigmadEprod = 0.; << 310 {//Probably that here the Bragg Model should be also used for kinEnergyProj/nuc<2MeV
286 G4double Emax_proj = GetSecondAdjEnergyMa << 311
287 G4double Emin_proj = GetSecondAdjEnergyMi << 312
288 << 313
289 // the produced particle should have a kinet << 314 G4double dSigmadEprod=0;
290 // projectile << 315 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
291 if(kinEnergyProj > Emin_proj && kinEnergyPro << 316 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
292 { << 317
293 G4double Tmax = kinEnergyProj; << 318
294 G4double E1 = kinEnergyProd; << 319 if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ //the produced particle should have a kinetic energy smaller than the projectile
295 //1.0006 factor seems to give the best dif << 320 G4double Tmax=kinEnergyProj;
296 G4double E2 = kinEnergyProd *1.0006; << 321
297 G4double sigma1, sigma2; << 322 G4double E1=kinEnergyProd;
298 if(kinEnergyProj > 2. * MeV) << 323 G4double E2=kinEnergyProd*1.000001;
299 { << 324 G4double dE=(E2-E1);
300 sigma1 = fDirectModel->ComputeCrossSecti << 325 G4double sigma1,sigma2;
301 fDirectPrimaryPart, kinEnergyProj, Z, << 326 if (kinEnergyProj >2.*MeV){
302 sigma2 = fDirectModel->ComputeCrossSecti << 327 sigma1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
303 fDirectPrimaryPart, kinEnergyProj, Z, << 328 sigma2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
304 } << 329 }
305 else << 330 else {
306 { << 331 sigma1=theBraggDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
307 sigma1 = fBraggDirectEMModel->ComputeCro << 332 sigma2=theBraggDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
308 fDirectPrimaryPart, kinEnergyProj, Z, << 333 }
309 sigma2 = fBraggDirectEMModel->ComputeCro << 334
310 fDirectPrimaryPart, kinEnergyProj, Z, << 335
311 } << 336 dSigmadEprod=(sigma1-sigma2)/dE;
312 << 337 if (dSigmadEprod>1.) {
313 dSigmadEprod = (sigma1 - sigma2) / (E2 - E << 338 G4cout<<"sigma1 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma1<<G4endl;
314 if(dSigmadEprod > 1.) << 339 G4cout<<"sigma2 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma2<<G4endl;
315 { << 340 G4cout<<"dsigma "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<dSigmadEprod<<G4endl;
316 G4cout << "sigma1 " << kinEnergyProj / M << 341
317 << '\t' << sigma1 << G4endl; << 342 }
318 G4cout << "sigma2 " << kinEnergyProj / M << 343
319 << '\t' << sigma2 << G4endl; << 344
320 G4cout << "dsigma " << kinEnergyProj / M << 345
321 << '\t' << dSigmadEprod << G4endl << 346 //correction of differential cross section at high energy to correct for the suppression of particle at secondary at high
322 } << 347 //energy used in the Bethe Bloch Model. This correction consist to multiply by g the probability function used
323 << 348 //to test the rejection of a secondary
324 // correction of differential cross sectio << 349 //-------------------------
325 // the suppression of particle at secondar << 350
326 // Bloch Model. This correction consists o << 351 //Source code taken from G4BetheBlochModel::SampleSecondaries
327 // function used to test the rejection of << 352
328 // from G4BetheBlochModel::SampleSecondari << 353 G4double deltaKinEnergy = kinEnergyProd;
329 G4double deltaKinEnergy = kinEnergyProd; << 354
330 << 355 //Part of the taken code
331 // projectile formfactor - suppression of << 356 //----------------------
332 // delta-electron production at high energ << 357
333 G4double x = fFormFact * deltaKinEnergy; << 358
334 if(x > 1.e-6) << 359
335 { << 360 // projectile formfactor - suppresion of high energy
336 G4double totEnergy = kinEnergyProj + fMa << 361 // delta-electron production at high energy
337 G4double etot2 = totEnergy * totEner << 362 G4double x = formfact*deltaKinEnergy;
338 G4double beta2 = kinEnergyProj * (kinEne << 363 if(x > 1.e-6) {
339 G4double f = 1.0 - beta2 * deltaKinE << 364
340 G4double f1 = 0.0; << 365
341 if(0.5 == fSpin) << 366 G4double totEnergy = kinEnergyProj + mass;
342 { << 367 G4double etot2 = totEnergy*totEnergy;
343 f1 = 0.5 * deltaKinEnergy * deltaKinEn << 368 G4double beta2 = kinEnergyProj*(kinEnergyProj + 2.0*mass)/etot2;
344 f += f1; << 369 G4double f;
345 } << 370 G4double f1 = 0.0;
346 G4double x1 = 1.0 + x; << 371 f = 1.0 - beta2*deltaKinEnergy/Tmax;
347 G4double gg = 1.0 / (x1 * x1); << 372 if( 0.5 == spin ) {
348 if(0.5 == fSpin) << 373 f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
349 { << 374 f += f1;
350 G4double x2 = 0.5 * electron_mass_c2 * << 375 }
351 gg *= (1.0 + fMagMoment2 * (x2 - f1 / << 376 G4double x1 = 1.0 + x;
352 } << 377 G4double gg = 1.0/(x1*x1);
353 if(gg > 1.0) << 378 if( 0.5 == spin ) {
354 { << 379 G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
355 G4cout << "### G4BetheBlochModel in Ad << 380 gg *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
356 << G4endl; << 381 }
357 gg = 1.; << 382 if(gg > 1.0) {
358 } << 383 G4cout << "### G4BetheBlochModel in Adjoint Sim WARNING: g= " << g
359 dSigmadEprod *= gg; << 384 << G4endl;
360 } << 385 gg=1.;
>> 386 }
>> 387 //G4cout<<"gg"<<gg<<G4endl;
>> 388 dSigmadEprod*=gg;
>> 389 }
>> 390
361 } 391 }
362 392
363 return dSigmadEprod; << 393 return dSigmadEprod;
364 } 394 }
365 395
366 ////////////////////////////////////////////// << 396
>> 397
>> 398 //////////////////////////////////////////////////////////////////////////////////////////////
>> 399 //
367 void G4AdjointhIonisationModel::DefineProjecti 400 void G4AdjointhIonisationModel::DefineProjectileProperty()
368 { << 401 {
369 // Slightly modified code taken from G4Bethe << 402 //Slightly modified code taken from G4BetheBlochModel::SetParticle
370 G4String pname = fDirectPrimaryPart->GetPart << 403 //------------------------------------------------
371 << 404 G4String pname = theDirectPrimaryPartDef->GetParticleName();
372 fMass = fDirectPrimaryPart->GetPDG << 405 if (theDirectPrimaryPartDef->GetParticleType() == "nucleus" &&
373 fSpin = fDirectPrimaryPart->GetPDG << 406 pname != "deuteron" && pname != "triton") {
374 fMassRatio = electron_mass_c2 / fMass; << 407 isIon = true;
375 fOnePlusRatio2 = (1. + fMassRatio) * (1. + <<
376 fOneMinusRatio2 = (1. - fMassRatio) * (1. - <<
377 G4double magmom = fDirectPrimaryPart->GetPDG <<
378 (0.5 * eplus * hbar_Planck <<
379 fMagMoment2 = magmom * magmom - 1.0; <<
380 fFormFact = 0.0; <<
381 if(fDirectPrimaryPart->GetLeptonNumber() == <<
382 { <<
383 G4double x = 0.8426 * GeV; <<
384 if(fSpin == 0.0 && fMass < GeV) <<
385 { <<
386 x = 0.736 * GeV; <<
387 } <<
388 else if(fMass > GeV) <<
389 { <<
390 x /= G4NistManager::Instance()->GetZ13(f <<
391 } 408 }
392 fFormFact = 2.0 * electron_mass_c2 / (x * << 409
393 } << 410 mass = theDirectPrimaryPartDef->GetPDGMass();
>> 411 mass_ratio_projectile = proton_mass_c2/theDirectPrimaryPartDef->GetPDGMass();;
>> 412 spin = theDirectPrimaryPartDef->GetPDGSpin();
>> 413 G4double q = theDirectPrimaryPartDef->GetPDGCharge()/eplus;
>> 414 chargeSquare = q*q;
>> 415 ratio = electron_mass_c2/mass;
>> 416 ratio2 = ratio*ratio;
>> 417 one_plus_ratio_2=(1+ratio)*(1+ratio);
>> 418 one_minus_ratio_2=(1-ratio)*(1-ratio);
>> 419 G4double magmom = theDirectPrimaryPartDef->GetPDGMagneticMoment()
>> 420 *mass/(0.5*eplus*hbar_Planck*c_squared);
>> 421 magMoment2 = magmom*magmom - 1.0;
>> 422 formfact = 0.0;
>> 423 if(theDirectPrimaryPartDef->GetLeptonNumber() == 0) {
>> 424 G4double x = 0.8426*GeV;
>> 425 if(spin == 0.0 && mass < GeV) {x = 0.736*GeV;}
>> 426 else if(mass > GeV) {
>> 427 x /= G4NistManager::Instance()->GetZ13(mass/proton_mass_c2);
>> 428 // tlimit = 51.2*GeV*A13[iz]*A13[iz];
>> 429 }
>> 430 formfact = 2.0*electron_mass_c2/(x*x);
>> 431 tlimit = 2.0/formfact;
>> 432 }
394 } 433 }
395 434
396 ////////////////////////////////////////////// 435 ////////////////////////////////////////////////////////////////////////////////
397 G4double G4AdjointhIonisationModel::AdjointCro << 436 //
398 const G4MaterialCutsCouple* aCouple, G4doubl << 437 G4double G4AdjointhIonisationModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
399 G4bool isScatProjToProj) << 438 G4double primEnergy,
400 { << 439 G4bool IsScatProjToProjCase)
401 if(fUseMatrix) << 440 {
402 return G4VEmAdjointModel::AdjointCrossSect << 441 if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase);
403 <<
404 DefineCurrentMaterial(aCouple); 442 DefineCurrentMaterial(aCouple);
>> 443
>> 444
>> 445 G4double Cross=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2*mass;
>> 446
>> 447 if (!IsScatProjToProjCase ){
>> 448 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy);
>> 449 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy);
>> 450 if (Emax_proj>Emin_proj && primEnergy > currentTcutForDirectSecond) {
>> 451 Cross*=(1./Emin_proj -1./Emax_proj)/primEnergy;
>> 452 }
>> 453 else Cross=0.;
>> 454
>> 455
>> 456
>> 457
>> 458
>> 459
>> 460 }
>> 461 else {
>> 462 G4double Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy);
>> 463 G4double Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,currentTcutForDirectSecond);
>> 464 G4double diff1=Emin_proj-primEnergy;
>> 465 G4double diff2=Emax_proj-primEnergy;
>> 466 G4double t1=(1./diff1+1./Emin_proj-1./diff2-1./Emax_proj)/primEnergy;
>> 467 //G4double t2=2.*std::log(diff2*Emin_proj/Emax_proj/diff1)/primEnergy/primEnergy;
>> 468 G4double t2=2.*std::log(Emax_proj/Emin_proj)/primEnergy/primEnergy;
>> 469 Cross*=(t1+t2);
405 470
406 G4double Cross = << 471
407 fCurrentMaterial->GetElectronDensity() * t <<
408 <<
409 if(!isScatProjToProj) <<
410 { <<
411 G4double Emax_proj = GetSecondAdjEnergyMax <<
412 G4double Emin_proj = GetSecondAdjEnergyMin <<
413 if(Emax_proj > Emin_proj && primEnergy > f <<
414 { <<
415 Cross *= (1. / Emin_proj - 1. / Emax_pro <<
416 } <<
417 else <<
418 Cross = 0.; <<
419 } <<
420 else <<
421 { <<
422 G4double Emax_proj = GetSecondAdjEnergyMax <<
423 G4double Emin_proj = <<
424 GetSecondAdjEnergyMinForScatProjToProj(p <<
425 G4double diff1 = Emin_proj - primEnergy; <<
426 G4double diff2 = Emax_proj - primEnergy; <<
427 G4double t1 = <<
428 (1. / diff1 + 1. / Emin_proj - 1. / diff <<
429 G4double t2 = <<
430 2. * std::log(Emax_proj / Emin_proj) / p <<
431 Cross *= (t1 + t2); <<
432 } 472 }
433 fLastCS = Cross; << 473 lastCS =Cross;
434 return Cross; << 474 return Cross;
435 } 475 }
436 <<
437 ////////////////////////////////////////////// 476 //////////////////////////////////////////////////////////////////////////////
438 G4double G4AdjointhIonisationModel::GetSecondA << 477 //
439 G4double primAdjEnergy) << 478 G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
440 { << 479 {
441 G4double Tmax = primAdjEnergy * fOnePlusRati << 480 G4double Tmax=PrimAdjEnergy*one_plus_ratio_2/(one_minus_ratio_2-2.*ratio*PrimAdjEnergy/mass);
442 (fOneMinusRatio2 - 2. * fMas <<
443 return Tmax; 481 return Tmax;
444 } 482 }
445 <<
446 ////////////////////////////////////////////// 483 //////////////////////////////////////////////////////////////////////////////
447 G4double G4AdjointhIonisationModel::GetSecondA << 484 //
448 G4double primAdjEnergy, G4double tcut) << 485 G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy,G4double Tcut)
449 { << 486 { return PrimAdjEnergy+Tcut;
450 return primAdjEnergy + tcut; <<
451 } 487 }
452 <<
453 ////////////////////////////////////////////// 488 //////////////////////////////////////////////////////////////////////////////
454 G4double G4AdjointhIonisationModel::GetSecondA << 489 //
455 { << 490 G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMaxForProdToProjCase(G4double )
456 return GetHighEnergyLimit(); << 491 { return HighEnergyLimit;
457 } 492 }
458 <<
459 ////////////////////////////////////////////// 493 //////////////////////////////////////////////////////////////////////////////
460 G4double G4AdjointhIonisationModel::GetSecondA << 494 //
461 G4double primAdjEnergy) << 495 G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
462 { << 496 { G4double Tmin= (2*PrimAdjEnergy-4*mass + std::sqrt(4.*PrimAdjEnergy*PrimAdjEnergy +16.*mass*mass + 8.*PrimAdjEnergy*mass*(1/ratio +ratio)))/4.;
463 G4double Tmin = << 497 return Tmin;
464 (2. * primAdjEnergy - 4. * fMass + <<
465 std::sqrt(4. * primAdjEnergy * primAdjEne <<
466 8. * primAdjEnergy * fMass * (1 <<
467 4.; <<
468 return Tmin; <<
469 } 498 }
470 499