source: trunk/source/processes/electromagnetic/adjoint/src/G4AdjointhIonisationModel.cc@ 1198

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

nvx fichiers dans CVS
File size: 18.0 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: G4AdjointhIonisationModel.cc,v 1.2 2009/11/20 10:31:20 ldesorgh Exp $
27// GEANT4 tag $Name: geant4-09-03-cand-01 $
28//
29#include "G4AdjointhIonisationModel.hh"
30#include "G4AdjointCSManager.hh"
31#include "G4Integrator.hh"
32#include "G4TrackStatus.hh"
33#include "G4ParticleChange.hh"
34#include "G4AdjointElectron.hh"
35#include "G4AdjointProton.hh"
36#include "G4AdjointInterpolator.hh"
37#include "G4BetheBlochModel.hh"
38#include "G4BraggModel.hh"
39#include "G4Proton.hh"
40#include "G4NistManager.hh"
41
42////////////////////////////////////////////////////////////////////////////////
43//
44G4AdjointhIonisationModel::G4AdjointhIonisationModel(G4ParticleDefinition* projectileDefinition):
45 G4VEmAdjointModel("Adjoint_hIonisation")
46{
47
48
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
71 DefineProjectileProperty();
72
73
74
75
76
77
78
79
80}
81////////////////////////////////////////////////////////////////////////////////
82//
83G4AdjointhIonisationModel::~G4AdjointhIonisationModel()
84{;}
85
86
87////////////////////////////////////////////////////////////////////////////////
88//
89void G4AdjointhIonisationModel::SampleSecondaries(const G4Track& aTrack,
90 G4bool IsScatProjToProjCase,
91 G4ParticleChange* fParticleChange)
92{
93 if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange);
94
95 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
96
97 //Elastic inverse scattering
98 //---------------------------------------------------------
99 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
100 G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
101
102 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
103 return;
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
160
161}
162
163////////////////////////////////////////////////////////////////////////////////
164//
165void G4AdjointhIonisationModel::RapidSampleSecondaries(const G4Track& aTrack,
166 G4bool IsScatProjToProjCase,
167 G4ParticleChange* fParticleChange)
168{
169
170 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
171 DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
172
173
174 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
175 G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
176
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.*log(f32);
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=adjointPrimKinEnergy/(1.-f31*pow(f32,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
301
302
303
304
305
306
307}
308
309////////////////////////////////////////////////////////////////////////////////
310//
311G4double G4AdjointhIonisationModel::DiffCrossSectionPerAtomPrimToSecond(
312 G4double kinEnergyProj,
313 G4double kinEnergyProd,
314 G4double Z,
315 G4double A)
316{//Probably that here the Bragg Model should be also used for kinEnergyProj/nuc<2MeV
317
318
319
320 G4double dSigmadEprod=0;
321 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
322 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
323
324
325 if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ //the produced particle should have a kinetic energy smaller than the projectile
326 G4double Tmax=kinEnergyProj;
327
328 G4double E1=kinEnergyProd;
329 G4double E2=kinEnergyProd*1.000001;
330 G4double dE=(E2-E1);
331 G4double sigma1,sigma2;
332 if (kinEnergyProj >2.*MeV){
333 sigma1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
334 sigma2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
335 }
336 else {
337 sigma1=theBraggDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
338 sigma2=theBraggDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
339 }
340
341
342 dSigmadEprod=(sigma1-sigma2)/dE;
343 if (dSigmadEprod>1.) {
344 G4cout<<"sigma1 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma1<<G4endl;
345 G4cout<<"sigma2 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma2<<G4endl;
346 G4cout<<"dsigma "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<dSigmadEprod<<G4endl;
347
348 }
349
350
351
352 //correction of differential cross section at high energy to correct for the suppression of particle at secondary at high
353 //energy used in the Bethe Bloch Model. This correction consist to multiply by g the probability function used
354 //to test the rejection of a secondary
355 //-------------------------
356
357 //Source code taken from G4BetheBlochModel::SampleSecondaries
358
359 G4double deltaKinEnergy = kinEnergyProd;
360
361 //Part of the taken code
362 //----------------------
363
364
365
366 // projectile formfactor - suppresion of high energy
367 // delta-electron production at high energy
368 G4double x = formfact*deltaKinEnergy;
369 if(x > 1.e-6) {
370
371
372 G4double totEnergy = kinEnergyProj + mass;
373 G4double etot2 = totEnergy*totEnergy;
374 G4double beta2 = kinEnergyProj*(kinEnergyProj + 2.0*mass)/etot2;
375 G4double f;
376 G4double f1 = 0.0;
377 f = 1.0 - beta2*deltaKinEnergy/Tmax;
378 if( 0.5 == spin ) {
379 f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
380 f += f1;
381 }
382 G4double x1 = 1.0 + x;
383 G4double g = 1.0/(x1*x1);
384 if( 0.5 == spin ) {
385 G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
386 g *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
387 }
388 if(g > 1.0) {
389 G4cout << "### G4BetheBlochModel in Adjoint Sim WARNING: g= " << g
390 << G4endl;
391 g=1.;
392 }
393 //G4cout<<"g"<<g<<G4endl;
394 dSigmadEprod*=g;
395 }
396
397 }
398
399 return dSigmadEprod;
400}
401
402
403
404//////////////////////////////////////////////////////////////////////////////////////////////
405//
406void G4AdjointhIonisationModel::DefineProjectileProperty()
407{
408 //Slightly modified code taken from G4BetheBlochModel::SetParticle
409 //------------------------------------------------
410 G4String pname = theDirectPrimaryPartDef->GetParticleName();
411 if (theDirectPrimaryPartDef->GetParticleType() == "nucleus" &&
412 pname != "deuteron" && pname != "triton") {
413 isIon = true;
414 }
415
416 mass = theDirectPrimaryPartDef->GetPDGMass();
417 spin = theDirectPrimaryPartDef->GetPDGSpin();
418 G4double q = theDirectPrimaryPartDef->GetPDGCharge()/eplus;
419 chargeSquare = q*q;
420 ratio = electron_mass_c2/mass;
421 ratio2 = ratio*ratio;
422 one_plus_ratio_2=(1+ratio)*(1+ratio);
423 one_minus_ratio_2=(1-ratio)*(1-ratio);
424 G4double magmom = theDirectPrimaryPartDef->GetPDGMagneticMoment()
425 *mass/(0.5*eplus*hbar_Planck*c_squared);
426 magMoment2 = magmom*magmom - 1.0;
427 formfact = 0.0;
428 if(theDirectPrimaryPartDef->GetLeptonNumber() == 0) {
429 G4double x = 0.8426*GeV;
430 if(spin == 0.0 && mass < GeV) {x = 0.736*GeV;}
431 else if(mass > GeV) {
432 x /= G4NistManager::Instance()->GetZ13(mass/proton_mass_c2);
433 // tlimit = 51.2*GeV*A13[iz]*A13[iz];
434 }
435 formfact = 2.0*electron_mass_c2/(x*x);
436 tlimit = 2.0/formfact;
437 }
438}
439
440////////////////////////////////////////////////////////////////////////////////
441//
442G4double G4AdjointhIonisationModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
443 G4double primEnergy,
444 G4bool IsScatProjToProjCase)
445{
446 if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase);
447 DefineCurrentMaterial(aCouple);
448
449
450 G4double Cross=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2*mass;
451
452 if (!IsScatProjToProjCase ){
453 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy);
454 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy);
455 if (Emax_proj>Emin_proj && primEnergy > currentTcutForDirectSecond) {
456 Cross*=(1./Emin_proj -1./Emax_proj)/primEnergy;
457 }
458 else Cross=0.;
459
460
461
462
463
464
465 }
466 else {
467 G4double Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy);
468 G4double Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,currentTcutForDirectSecond);
469 G4double diff1=Emin_proj-primEnergy;
470 G4double diff2=Emax_proj-primEnergy;
471 G4double t1=(1./diff1+1./Emin_proj-1./diff2-1./Emax_proj)/primEnergy;
472 G4double t2=2.*log(diff2*Emin_proj/Emax_proj/diff1)/primEnergy/primEnergy;
473 Cross*=(t1+t2);
474
475
476 }
477 lastCS =Cross;
478 return Cross;
479}
480//////////////////////////////////////////////////////////////////////////////
481//
482G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
483{
484 G4double Tmax=PrimAdjEnergy*one_plus_ratio_2/(one_minus_ratio_2-2.*ratio*PrimAdjEnergy/mass);
485 return Tmax;
486}
487//////////////////////////////////////////////////////////////////////////////
488//
489G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy,G4double Tcut)
490{ return PrimAdjEnergy+Tcut;
491}
492//////////////////////////////////////////////////////////////////////////////
493//
494G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMaxForProdToProjCase(G4double )
495{ return HighEnergyLimit;
496}
497//////////////////////////////////////////////////////////////////////////////
498//
499G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
500{ G4double Tmin= (2*PrimAdjEnergy-4*mass + std::sqrt(4.*PrimAdjEnergy*PrimAdjEnergy +16.*mass*mass + 8.*PrimAdjEnergy*mass*(1/ratio +ratio)))/4.;
501 return Tmin;
502}
Note: See TracBrowser for help on using the repository browser.