source: trunk/source/processes/electromagnetic/adjoint/src/G4AdjointIonIonisationModel.cc@ 1201

Last change on this file since 1201 was 1197, checked in by garnier, 16 years ago

nvx fichiers dans CVS
File size: 14.7 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// $Id: G4AdjointIonIonisationModel.cc,v 1.2 2009/11/20 10:31:20 ldesorgh Exp $
27// GEANT4 tag $Name: geant4-09-03-cand-01 $
28//
29#include "G4AdjointIonIonisationModel.hh"
30#include "G4AdjointCSManager.hh"
31
32
33#include "G4Integrator.hh"
34#include "G4TrackStatus.hh"
35#include "G4ParticleChange.hh"
36#include "G4AdjointElectron.hh"
37#include "G4AdjointProton.hh"
38#include "G4AdjointInterpolator.hh"
39#include "G4BetheBlochModel.hh"
40#include "G4BraggIonModel.hh"
41#include "G4Proton.hh"
42#include "G4GenericIon.hh"
43#include "G4NistManager.hh"
44
45////////////////////////////////////////////////////////////////////////////////
46//
47G4AdjointIonIonisationModel::G4AdjointIonIonisationModel():
48 G4VEmAdjointModel("Adjoint_IonIonisation")
49{
50
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 use_only_bragg = false; // for the Ion ionisation using the parametrised table model the cross sections and the sample of secondaries is done
59 // as in the BraggIonModel, Therefore the use of this flag;
60
61 //The direct EM Model is taken has BetheBloch it is only used for the computation
62 // of the differential cross section.
63 //The Bragg model could be used as an alternative as it offers the same differential cross section
64
65 theBetheBlochDirectEMModel = new G4BetheBlochModel(G4GenericIon::GenericIon());
66 theBraggIonDirectEMModel = new G4BraggIonModel(G4GenericIon::GenericIon());
67 theAdjEquivOfDirectSecondPartDef=G4AdjointElectron::AdjointElectron();
68 theDirectPrimaryPartDef =0;
69 theAdjEquivOfDirectPrimPartDef =0;
70 /* theDirectPrimaryPartDef =fwd_ion;
71 theAdjEquivOfDirectPrimPartDef =adj_ion;
72
73 DefineProjectileProperty();*/
74
75
76
77
78}
79////////////////////////////////////////////////////////////////////////////////
80//
81G4AdjointIonIonisationModel::~G4AdjointIonIonisationModel()
82{;}
83////////////////////////////////////////////////////////////////////////////////
84//
85void G4AdjointIonIonisationModel::SampleSecondaries(const G4Track& aTrack,
86 G4bool IsScatProjToProjCase,
87 G4ParticleChange* fParticleChange)
88{
89 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
90
91 //Elastic inverse scattering
92 //---------------------------------------------------------
93 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
94 G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
95
96 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
97 return;
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
154
155}
156
157////////////////////////////////////////////////////////////////////////////////
158//
159G4double G4AdjointIonIonisationModel::DiffCrossSectionPerAtomPrimToSecond(
160 G4double kinEnergyProj,
161 G4double kinEnergyProd,
162 G4double Z,
163 G4double A)
164{//Probably that here the Bragg Model should be also used for kinEnergyProj/nuc<2MeV
165
166
167
168 G4double dSigmadEprod=0;
169 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
170 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
171
172 G4double kinEnergyProjScaled = massRatio*kinEnergyProj;
173
174
175 if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ //the produced particle should have a kinetic energy smaller than the projectile
176 G4double Tmax=kinEnergyProj;
177
178 G4double E1=kinEnergyProd;
179 G4double E2=kinEnergyProd*1.000001;
180 G4double dE=(E2-E1);
181 G4double sigma1,sigma2;
182 theDirectEMModel =theBraggIonDirectEMModel;
183 if (kinEnergyProjScaled >2.*MeV && !use_only_bragg) theDirectEMModel = theBetheBlochDirectEMModel; //Bethe Bloch Model
184 sigma1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
185 sigma2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
186
187 dSigmadEprod=(sigma1-sigma2)/dE;
188
189 //G4double chargeSqRatio = currentModel->GetChargeSquareRatio(theDirectPrimaryPartDef,currentMaterial,E);
190
191
192
193 if (dSigmadEprod>1.) {
194 G4cout<<"sigma1 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma1<<G4endl;
195 G4cout<<"sigma2 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma2<<G4endl;
196 G4cout<<"dsigma "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<dSigmadEprod<<G4endl;
197
198 }
199
200
201
202
203
204
205 if (theDirectEMModel == theBetheBlochDirectEMModel ){
206 //correction of differential cross section at high energy to correct for the suppression of particle at secondary at high
207 //energy used in the Bethe Bloch Model. This correction consist to multiply by g the probability function used
208 //to test the rejection of a secondary
209 //-------------------------
210
211 //Source code taken from G4BetheBlochModel::SampleSecondaries
212
213 G4double deltaKinEnergy = kinEnergyProd;
214
215 //Part of the taken code
216 //----------------------
217
218
219
220 // projectile formfactor - suppresion of high energy
221 // delta-electron production at high energy
222
223
224 G4double x = formfact*deltaKinEnergy;
225 if(x > 1.e-6) {
226 G4double totEnergy = kinEnergyProj + mass;
227 G4double etot2 = totEnergy*totEnergy;
228 G4double beta2 = kinEnergyProj*(kinEnergyProj + 2.0*mass)/etot2;
229 G4double f;
230 G4double f1 = 0.0;
231 f = 1.0 - beta2*deltaKinEnergy/Tmax;
232 if( 0.5 == spin ) {
233 f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
234 f += f1;
235 }
236 G4double x1 = 1.0 + x;
237 G4double g = 1.0/(x1*x1);
238 if( 0.5 == spin ) {
239 G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
240 g *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
241 }
242 if(g > 1.0) {
243 G4cout << "### G4BetheBlochModel in Adjoint Sim WARNING: g= " << g
244 << G4endl;
245 g=1.;
246 }
247 //G4cout<<"g"<<g<<G4endl;
248 dSigmadEprod*=g;
249 }
250 }
251
252 }
253
254 return dSigmadEprod;
255}
256//////////////////////////////////////////////////////////////////////////////////////////////
257//
258void G4AdjointIonIonisationModel::SetIon(G4ParticleDefinition* adj_ion, G4ParticleDefinition* fwd_ion)
259{ theDirectPrimaryPartDef =fwd_ion;
260 theAdjEquivOfDirectPrimPartDef =adj_ion;
261
262 DefineProjectileProperty();
263}
264//////////////////////////////////////////////////////////////////////////////////////////////
265//
266void G4AdjointIonIonisationModel::CorrectPostStepWeight(G4ParticleChange* fParticleChange, G4double old_weight,
267 G4double adjointPrimKinEnergy, G4double projectileKinEnergy,G4bool )
268{
269 //It is needed because the direct cross section used to compute the differential cross section is not the one used in
270 // the direct model where the GenericIon stuff is considered with correction of effective charge. In the direct model the samnepl of secondaries does
271 // not reflect the integral cross section. The integral fwd cross section that we used to compute the differential CS
272 // match the sample of secondaries in the forward case despite the fact that its is not the same total CS than in the FWD case. For this reasion an extra
273 // weight correction is needed at the end.
274
275
276 G4double new_weight=old_weight;
277
278 //the correction of CS due to the problem explained above
279 G4double kinEnergyProjScaled = massRatio*projectileKinEnergy;
280 theDirectEMModel =theBraggIonDirectEMModel;
281 if (kinEnergyProjScaled >2.*MeV && !use_only_bragg) theDirectEMModel = theBetheBlochDirectEMModel; //Bethe Bloch Model
282 G4double UsedFwdCS=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,projectileKinEnergy,1,1 ,currentTcutForDirectSecond,1.e20);
283 G4double chargeSqRatio =1.;
284 if (chargeSquare>1.) chargeSqRatio = theDirectEMModel->GetChargeSquareRatio(theDirectPrimaryPartDef,currentMaterial,projectileKinEnergy);
285 G4double CorrectFwdCS = chargeSqRatio*theDirectEMModel->ComputeCrossSectionPerAtom(G4GenericIon::GenericIon(),kinEnergyProjScaled,1,1 ,currentTcutForDirectSecond,1.e20);
286 if (UsedFwdCS >0) new_weight*= CorrectFwdCS/UsedFwdCS;//May be some check is needed if UsedFwdCS ==0 probably that then we should avoid a secondary to be produced,
287
288
289 //additional CS crorrection needed for cross section biasing in general.
290 //May be wrong for ions!!! Most of the time not used!
291 G4double w_corr =1./CS_biasing_factor;
292 w_corr*=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
293 new_weight*=w_corr;
294
295 new_weight*=projectileKinEnergy/adjointPrimKinEnergy;
296
297 fParticleChange->SetParentWeightByProcess(false);
298 fParticleChange->SetSecondaryWeightByProcess(false);
299 fParticleChange->ProposeParentWeight(new_weight);
300}
301
302
303//////////////////////////////////////////////////////////////////////////////////////////////
304//
305void G4AdjointIonIonisationModel::DefineProjectileProperty()
306{
307 //Slightly modified code taken from G4BetheBlochModel::SetParticle
308 //------------------------------------------------
309 G4String pname = theDirectPrimaryPartDef->GetParticleName();
310 if (theDirectPrimaryPartDef->GetParticleType() == "nucleus" &&
311 pname != "deuteron" && pname != "triton") {
312 isIon = true;
313 }
314
315 mass = theDirectPrimaryPartDef->GetPDGMass();
316 massRatio= G4GenericIon::GenericIon()->GetPDGMass()/mass;
317 spin = theDirectPrimaryPartDef->GetPDGSpin();
318 G4double q = theDirectPrimaryPartDef->GetPDGCharge()/eplus;
319 chargeSquare = q*q;
320 ratio = electron_mass_c2/mass;
321 ratio2 = ratio*ratio;
322 one_plus_ratio_2=(1+ratio)*(1+ratio);
323 one_minus_ratio_2=(1-ratio)*(1-ratio);
324 G4double magmom = theDirectPrimaryPartDef->GetPDGMagneticMoment()
325 *mass/(0.5*eplus*hbar_Planck*c_squared);
326 magMoment2 = magmom*magmom - 1.0;
327 formfact = 0.0;
328 if(theDirectPrimaryPartDef->GetLeptonNumber() == 0) {
329 G4double x = 0.8426*GeV;
330 if(spin == 0.0 && mass < GeV) {x = 0.736*GeV;}
331 else if(mass > GeV) {
332 x /= G4NistManager::Instance()->GetZ13(mass/proton_mass_c2);
333 // tlimit = 51.2*GeV*A13[iz]*A13[iz];
334 }
335 formfact = 2.0*electron_mass_c2/(x*x);
336 tlimit = 2.0/formfact;
337 }
338}
339
340
341//////////////////////////////////////////////////////////////////////////////
342//
343G4double G4AdjointIonIonisationModel::GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
344{
345 G4double Tmax=PrimAdjEnergy*one_plus_ratio_2/(one_minus_ratio_2-2.*ratio*PrimAdjEnergy/mass);
346 return Tmax;
347}
348//////////////////////////////////////////////////////////////////////////////
349//
350G4double G4AdjointIonIonisationModel::GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy,G4double Tcut)
351{ return PrimAdjEnergy+Tcut;
352}
353//////////////////////////////////////////////////////////////////////////////
354//
355G4double G4AdjointIonIonisationModel::GetSecondAdjEnergyMaxForProdToProjCase(G4double )
356{ return HighEnergyLimit;
357}
358//////////////////////////////////////////////////////////////////////////////
359//
360G4double G4AdjointIonIonisationModel::GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
361{ G4double Tmin= (2*PrimAdjEnergy-4*mass + std::sqrt(4.*PrimAdjEnergy*PrimAdjEnergy +16.*mass*mass + 8.*PrimAdjEnergy*mass*(1/ratio +ratio)))/4.;
362 return Tmin;
363}
Note: See TracBrowser for help on using the repository browser.