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