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