source: trunk/source/processes/electromagnetic/adjoint/src/G4AdjointBremsstrahlungModel.cc @ 1315

Last change on this file since 1315 was 1228, checked in by garnier, 15 years ago

update geant4.9.3 tag
File size: 16.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: G4AdjointBremsstrahlungModel.cc,v 1.5 2009/12/16 17:50:01 gunter Exp $
27// GEANT4 tag $Name: geant4-09-03 $
28//
29#include "G4AdjointBremsstrahlungModel.hh"
30#include "G4AdjointCSManager.hh"
31#include "G4Integrator.hh"
32#include "G4TrackStatus.hh"
33#include "G4ParticleChange.hh"
34#include "G4AdjointElectron.hh"
35#include "G4AdjointGamma.hh"
36#include "G4Electron.hh"
37
38#include "G4Timer.hh"
39//#include "G4PenelopeBremsstrahlungModel.hh"
40
41
42////////////////////////////////////////////////////////////////////////////////
43//
44G4AdjointBremsstrahlungModel::G4AdjointBremsstrahlungModel():
45 G4VEmAdjointModel("AdjointeBremModel"),
46  MigdalConstant(classic_electr_radius*electron_Compton_length*electron_Compton_length*4.0*pi)
47{ 
48  SetUseMatrix(false);
49  SetUseMatrixPerElement(false);
50 
51  theDirectStdBremModel = new G4eBremsstrahlungModel(G4Electron::Electron(),"TheDirecteBremModel");
52  theDirectEMModel=theDirectStdBremModel;
53 // theDirectPenelopeBremModel =0;
54       
55  SetApplyCutInRange(true);
56  highKinEnergy= 100.*TeV;
57  lowKinEnergy = 1.0*keV;
58  theTimer =new G4Timer();
59 
60  theAdjEquivOfDirectPrimPartDef =G4AdjointElectron::AdjointElectron();
61  theAdjEquivOfDirectSecondPartDef=G4AdjointGamma::AdjointGamma();
62  theDirectPrimaryPartDef=G4Electron::Electron();
63  second_part_of_same_type=false;
64 
65  /*UsePenelopeModel=false;
66  if (UsePenelopeModel) {
67        G4PenelopeBremsstrahlungModel* thePenelopeModel = new G4PenelopeBremsstrahlungModel(G4Electron::Electron(),"PenelopeBrem");
68        theEmModelManagerForFwdModels = new G4EmModelManager();
69        isPenelopeModelInitialised = false;
70        G4VEmFluctuationModel* f=0;
71        G4Region* r=0;
72        theDirectEMModel=thePenelopeModel;
73        theEmModelManagerForFwdModels->AddEmModel(1, thePenelopeModel, f, r);
74  }
75  */   
76 
77
78 
79}
80
81////////////////////////////////////////////////////////////////////////////////
82//
83void G4AdjointBremsstrahlungModel::SampleSecondaries(const G4Track& aTrack,
84                       G4bool IsScatProjToProjCase,
85                       G4ParticleChange* fParticleChange)
86{
87 if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange); 
88
89 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
90 DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
91 
92 
93 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
94 G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy();
95 
96 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
97        return;
98 }
99 
100  G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy,
101                                                        IsScatProjToProjCase);
102 //Weight correction
103 //-----------------------                                         
104 CorrectPostStepWeight(fParticleChange, 
105                       aTrack.GetWeight(), 
106                       adjointPrimKinEnergy,
107                       projectileKinEnergy,
108                       IsScatProjToProjCase);   
109 
110 
111 //Kinematic
112 //---------
113 G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
114 G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
115 G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;       
116 G4double projectileP = std::sqrt(projectileP2);
117 
118 
119 //Angle of the gamma direction with the projectile taken from G4eBremsstrahlungModel
120 //------------------------------------------------
121  G4double u;
122  const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
123
124  if (9./(9.+d) > G4UniformRand()) u = - std::log(G4UniformRand()*G4UniformRand())/a1;
125     else                          u = - std::log(G4UniformRand()*G4UniformRand())/a2;
126
127  G4double theta = u*electron_mass_c2/projectileTotalEnergy;
128
129  G4double sint = std::sin(theta);
130  G4double cost = std::cos(theta);
131
132  G4double phi = twopi * G4UniformRand() ;
133 
134  G4ThreeVector projectileMomentum;
135  projectileMomentum=G4ThreeVector(std::cos(phi)*sint,std::sin(phi)*sint,cost)*projectileP; //gamma frame
136  if (IsScatProjToProjCase) {//the adjoint primary is the scattered e-
137        G4ThreeVector gammaMomentum = (projectileTotalEnergy-adjointPrimTotalEnergy)*G4ThreeVector(0.,0.,1.);
138        G4ThreeVector dirProd=projectileMomentum-gammaMomentum;
139        G4double cost1 = std::cos(dirProd.angle(projectileMomentum));
140        G4double sint1 =  std::sqrt(1.-cost1*cost1);
141        projectileMomentum=G4ThreeVector(std::cos(phi)*sint1,std::sin(phi)*sint1,cost1)*projectileP;
142 
143  }
144 
145  projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection());
146 
147 
148 
149  if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
150        fParticleChange->ProposeTrackStatus(fStopAndKill);
151        fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
152  }
153  else {
154        fParticleChange->ProposeEnergy(projectileKinEnergy);
155        fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
156       
157  }     
158} 
159////////////////////////////////////////////////////////////////////////////////
160//
161void G4AdjointBremsstrahlungModel::RapidSampleSecondaries(const G4Track& aTrack,
162                       G4bool IsScatProjToProjCase,
163                       G4ParticleChange* fParticleChange)
164{ 
165
166 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
167 DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
168 
169 
170 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
171 G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy();
172 
173 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
174        return;
175 }
176 
177 G4double projectileKinEnergy =0.;
178 G4double gammaEnergy=0.;
179 G4double diffCSUsed=0.; 
180 if (!IsScatProjToProjCase){
181        gammaEnergy=adjointPrimKinEnergy;
182        G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
183        G4double Emin=  GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);;
184        if (Emin>=Emax) return;
185        projectileKinEnergy=Emin*std::pow(Emax/Emin,G4UniformRand());
186        diffCSUsed=lastCZ/projectileKinEnergy;
187       
188 }
189 else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
190        G4double Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond);
191        if (Emin>=Emax) return;
192        G4double f1=(Emin-adjointPrimKinEnergy)/Emin;
193        G4double f2=(Emax-adjointPrimKinEnergy)/Emax/f1;
194        //G4cout<<"f1 and f2 "<<f1<<'\t'<<f2<<G4endl;
195        projectileKinEnergy=adjointPrimKinEnergy/(1.-f1*std::pow(f2,G4UniformRand()));
196        gammaEnergy=projectileKinEnergy-adjointPrimKinEnergy;
197        diffCSUsed=lastCZ*adjointPrimKinEnergy/projectileKinEnergy/gammaEnergy;
198       
199 }
200 
201 
202 
203                                                       
204 //Weight correction
205 //-----------------------
206 //First w_corr is set to the ratio between adjoint total CS and fwd total CS
207 G4double w_corr=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
208
209 //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
210 //Here we consider the true  diffCS as the one obtained by the numericla differentiation over Tcut of the direct CS, corrected by the Migdal term.
211 //Basically any other differential CS   diffCS could be used here (example Penelope).
212 
213 G4double diffCS = DiffCrossSectionPerVolumePrimToSecond(currentMaterial, projectileKinEnergy, gammaEnergy);
214 w_corr*=diffCS/diffCSUsed;
215           
216 G4double new_weight = aTrack.GetWeight()*w_corr;
217 fParticleChange->SetParentWeightByProcess(false);
218 fParticleChange->SetSecondaryWeightByProcess(false);
219 fParticleChange->ProposeParentWeight(new_weight);
220 
221 //Kinematic
222 //---------
223 G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
224 G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
225 G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;       
226 G4double projectileP = std::sqrt(projectileP2);
227 
228 
229 //Angle of the gamma direction with the projectile taken from G4eBremsstrahlungModel
230 //------------------------------------------------
231  G4double u;
232  const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
233
234  if (9./(9.+d) > G4UniformRand()) u = - std::log(G4UniformRand()*G4UniformRand())/a1;
235     else                          u = - std::log(G4UniformRand()*G4UniformRand())/a2;
236
237  G4double theta = u*electron_mass_c2/projectileTotalEnergy;
238
239  G4double sint = std::sin(theta);
240  G4double cost = std::cos(theta);
241
242  G4double phi = twopi * G4UniformRand() ;
243 
244  G4ThreeVector projectileMomentum;
245  projectileMomentum=G4ThreeVector(std::cos(phi)*sint,std::sin(phi)*sint,cost)*projectileP; //gamma frame
246  if (IsScatProjToProjCase) {//the adjoint primary is the scattered e-
247        G4ThreeVector gammaMomentum = (projectileTotalEnergy-adjointPrimTotalEnergy)*G4ThreeVector(0.,0.,1.);
248        G4ThreeVector dirProd=projectileMomentum-gammaMomentum;
249        G4double cost1 = std::cos(dirProd.angle(projectileMomentum));
250        G4double sint1 =  std::sqrt(1.-cost1*cost1);
251        projectileMomentum=G4ThreeVector(std::cos(phi)*sint1,std::sin(phi)*sint1,cost1)*projectileP;
252 
253  }
254 
255  projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection());
256 
257 
258 
259  if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
260        fParticleChange->ProposeTrackStatus(fStopAndKill);
261        fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
262  }
263  else {
264        fParticleChange->ProposeEnergy(projectileKinEnergy);
265        fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
266       
267  }     
268} 
269////////////////////////////////////////////////////////////////////////////////
270//
271G4AdjointBremsstrahlungModel::~G4AdjointBremsstrahlungModel()
272{;}
273
274////////////////////////////////////////////////////////////////////////////////
275//
276G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond(const G4Material* aMaterial,
277                                      G4double kinEnergyProj,  // kinetic energy of the primary particle before the interaction
278                                      G4double kinEnergyProd // kinetic energy of the secondary particle
279                                      )
280{/*if (UsePenelopeModel && !isPenelopeModelInitialised) {
281        theEmModelManagerForFwdModels->Initialise(G4Electron::Electron(),G4Gamma::Gamma(),1.,0);
282        isPenelopeModelInitialised =true;
283 }
284 */                                                             
285 return  DiffCrossSectionPerVolumePrimToSecondApproximated2(aMaterial,
286                                                                 kinEnergyProj, 
287                                                                 kinEnergyProd);
288 /*return G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToSecond(aMaterial,
289                                                                 kinEnergyProj,
290                                                                 kinEnergyProd);*/                                                                                                                               
291}                                     
292
293////////////////////////////////////////////////////////////////////////////////
294//
295G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecondApproximated1(
296                                      const G4Material* aMaterial,
297                                      G4double kinEnergyProj,  // kinetic energy of the primary particle before the interaction
298                                      G4double kinEnergyProd // kinetic energy of the secondary particle
299                                      )
300{
301 G4double dCrossEprod=0.;
302 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
303 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
304 
305 
306 //In this approximation we consider that the secondary gammas are sampled with 1/Egamma energy distribution
307 //This is what is applied in the discrete standard model before the  rejection test  that make a cooerction
308 //The application of the same rejection function is not possble here.
309 //The differentiation of the CS over Ecut does not produce neither a good differential CS. That is due to the
310 // fact that in the discrete model the differential CS and the integrated CS are both fitted but separatly and
311 // therefore do not allow a correct numerical differentiation of the integrated CS to get the differential one.
312 // In the future we plan to use the brem secondary spectra from the G4Penelope implementation
313 
314 if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){
315        G4double sigma=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,1.*keV);
316        dCrossEprod=sigma/kinEnergyProd/std::log(kinEnergyProj/keV);
317 }
318 return dCrossEprod;
319 
320}
321
322////////////////////////////////////////////////////////////////////////////////
323//
324G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecondApproximated2(
325                                      const G4Material* material,
326                                      G4double kinEnergyProj,  // kinetic energy of the primary particle before the interaction
327                                      G4double kinEnergyProd // kinetic energy of the secondary particle
328                                      )
329{
330 //In this approximation we derive the direct cross section over Tcut=gamma energy, en after apply the Migdla correction factor
331  //used in the direct model
332 
333 G4double dCrossEprod=0.;
334 
335 const G4ElementVector* theElementVector = material->GetElementVector();
336 const double* theAtomNumDensityVector = material->GetAtomicNumDensityVector();
337 G4double dum=0.;
338 G4double E1=kinEnergyProd,E2=kinEnergyProd*1.001;
339 G4double dE=E2-E1;
340 for (size_t i=0; i<material->GetNumberOfElements(); i++) { 
341        G4double C1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,(*theElementVector)[i]->GetZ(),dum ,E1);
342        G4double C2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,(*theElementVector)[i]->GetZ(),dum,E2);
343        dCrossEprod += theAtomNumDensityVector[i] * (C1-C2)/dE;
344   
345 }
346 
347 //Now the Migdal correction
348 
349 G4double totalEnergy = kinEnergyProj+electron_mass_c2 ;
350 G4double kp2 = MigdalConstant*totalEnergy*totalEnergy
351                                             *(material->GetElectronDensity());
352                                             
353 
354 G4double MigdalFactor = 1./(1.+kp2/(kinEnergyProd*kinEnergyProd)); // its seems that the factor used in the CS compuation i the direct
355                                                                    //model is different than the one used in the secondary sampling by a
356                                                                    //factor (1.+kp2) To be checked!
357 
358 dCrossEprod*=MigdalFactor;
359 return dCrossEprod;
360 
361}
362////////////////////////////////////////////////////////////////////////////////
363//
364G4double G4AdjointBremsstrahlungModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
365                                             G4double primEnergy,
366                                             G4bool IsScatProjToProjCase)
367{/* if (UsePenelopeModel && !isPenelopeModelInitialised) {
368        theEmModelManagerForFwdModels->Initialise(G4Electron::Electron(),G4Gamma::Gamma(),1.,0);
369        isPenelopeModelInitialised =true;
370  }
371  */
372  if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase);
373  DefineCurrentMaterial(aCouple);
374  G4double Cross=0.;
375  lastCZ=theDirectEMModel->CrossSectionPerVolume(aCouple->GetMaterial(),theDirectPrimaryPartDef,100.*MeV,100.*MeV/std::exp(1.));//this give the constant above
376 
377  if (!IsScatProjToProjCase ){
378        G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy);
379        G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy);
380        if (Emax_proj>Emin_proj && primEnergy > currentTcutForDirectSecond) Cross= lastCZ*std::log(Emax_proj/Emin_proj);
381  }
382  else {
383        G4double Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy);
384        G4double Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,currentTcutForDirectSecond);
385        if (Emax_proj>Emin_proj) Cross= lastCZ*std::log((Emax_proj-primEnergy)*Emin_proj/Emax_proj/(Emin_proj-primEnergy));
386       
387  }
388  return Cross; 
389}                                           
390
391
392
393
394
Note: See TracBrowser for help on using the repository browser.