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

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

update CVS release candidate geant4.9.3.01
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.4 2009/11/20 10:31:20 ldesorgh Exp $
27// GEANT4 tag $Name: geant4-09-03-cand-01 $
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 = - log(G4UniformRand()*G4UniformRand())/a1;
125 else u = - 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*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*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 = - log(G4UniformRand()*G4UniformRand())/a1;
235 else u = - 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/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*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*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.