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

Last change on this file since 1197 was 1197, checked in by garnier, 15 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.