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