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

Last change on this file since 1337 was 1337, checked in by garnier, 14 years ago

tag geant4.9.4 beta 1 + modifs locales
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-04-beta-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.