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