Geant4  10.03
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 100666 2016-10-31 10:27:00Z 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;
57  G4Region* r=0;
59 
60  SetApplyCutInRange(true);
61  highKinEnergy= 1.*GeV;
62  lowKinEnergy = 1.0*keV;
63 
64  lastCZ =0.;
65 
66 
71 
72 
74 
75 
76 }
78 //
80  G4VEmAdjointModel("AdjointeBremModel")
81 {
82  SetUseMatrix(false);
84 
90  G4Region* r=0;
92  // theDirectPenelopeBremModel =0;
93  SetApplyCutInRange(true);
94  highKinEnergy= 1.*GeV;
95  lowKinEnergy = 1.0*keV;
96  lastCZ =0.;
101 }
103 //
107 }
108 
110 //
112  G4bool IsScatProjToProjCase,
113  G4ParticleChange* fParticleChange)
114 {
115  if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange);
116 
117  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
119 
120 
121  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
122  G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy();
123 
124  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
125  return;
126  }
127 
128  G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy,
129  IsScatProjToProjCase);
130  //Weight correction
131  //-----------------------
132  CorrectPostStepWeight(fParticleChange,
133  aTrack.GetWeight(),
134  adjointPrimKinEnergy,
135  projectileKinEnergy,
136  IsScatProjToProjCase);
137 
138 
139  //Kinematic
140  //---------
142  G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
143  G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
144  G4double projectileP = std::sqrt(projectileP2);
145 
146 
147  //Angle of the gamma direction with the projectile taken from G4eBremsstrahlungModel
148  //------------------------------------------------
149  G4double u;
150  const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
151 
152  if (9./(9.+d) > G4UniformRand()) u = - std::log(G4UniformRand()*G4UniformRand())/a1;
153  else u = - std::log(G4UniformRand()*G4UniformRand())/a2;
154 
155  G4double theta = u*electron_mass_c2/projectileTotalEnergy;
156 
157  G4double sint = std::sin(theta);
158  G4double cost = std::cos(theta);
159 
160  G4double phi = twopi * G4UniformRand() ;
161 
162  G4ThreeVector projectileMomentum;
163  projectileMomentum=G4ThreeVector(std::cos(phi)*sint,std::sin(phi)*sint,cost)*projectileP; //gamma frame
164  if (IsScatProjToProjCase) {//the adjoint primary is the scattered e-
165  G4ThreeVector gammaMomentum = (projectileTotalEnergy-adjointPrimTotalEnergy)*G4ThreeVector(0.,0.,1.);
166  G4ThreeVector dirProd=projectileMomentum-gammaMomentum;
167  G4double cost1 = std::cos(dirProd.angle(projectileMomentum));
168  G4double sint1 = std::sqrt(1.-cost1*cost1);
169  projectileMomentum=G4ThreeVector(std::cos(phi)*sint1,std::sin(phi)*sint1,cost1)*projectileP;
170 
171  }
172 
173  projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection());
174 
175 
176 
177  if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
178  fParticleChange->ProposeTrackStatus(fStopAndKill);
179  fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
180  }
181  else {
182  fParticleChange->ProposeEnergy(projectileKinEnergy);
183  fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
184 
185  }
186 }
188 //
190  G4bool IsScatProjToProjCase,
191  G4ParticleChange* fParticleChange)
192 {
193 
194  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
196 
197 
198  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
199  G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy();
200 
201  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
202  return;
203  }
204 
205  G4double projectileKinEnergy =0.;
206  G4double gammaEnergy=0.;
207  G4double diffCSUsed=0.;
208  if (!IsScatProjToProjCase){
209  gammaEnergy=adjointPrimKinEnergy;
210  G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
211  G4double Emin= GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);;
212  if (Emin>=Emax) return;
213  projectileKinEnergy=Emin*std::pow(Emax/Emin,G4UniformRand());
214  diffCSUsed=CS_biasing_factor*lastCZ/projectileKinEnergy;
215 
216  }
217  else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
219  if (Emin>=Emax) return;
220  G4double f1=(Emin-adjointPrimKinEnergy)/Emin;
221  G4double f2=(Emax-adjointPrimKinEnergy)/Emax/f1;
222  projectileKinEnergy=adjointPrimKinEnergy/(1.-f1*std::pow(f2,G4UniformRand()));
223  gammaEnergy=projectileKinEnergy-adjointPrimKinEnergy;
224  diffCSUsed=lastCZ*adjointPrimKinEnergy/projectileKinEnergy/gammaEnergy;
225 
226  }
227 
228 
229 
230 
231  //Weight correction
232  //-----------------------
233  //First w_corr is set to the ratio between adjoint total CS and fwd total CS
234  //if this has to be done in the model
235  //For the case of forced interaction this will be done in the PostStepDoIt of the
236  //forced interaction
237  //It is important to set the weight before the vreation of the secondary
238  //
242  }
243  //G4cout<<"Correction factor start in brem model "<<w_corr<<std::endl;
244 
245 
246  //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
247  //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.
248  //Basically any other differential CS diffCS could be used here (example Penelope).
249 
250  G4double diffCS = DiffCrossSectionPerVolumePrimToSecond(currentMaterial, projectileKinEnergy, gammaEnergy);
251  /*G4cout<<"diffCS "<<diffCS <<std::endl;
252  G4cout<<"diffCS_Used "<<diffCSUsed <<std::endl;*/
253  w_corr*=diffCS/diffCSUsed;
254 
255 
256  G4double new_weight = aTrack.GetWeight()*w_corr;
257  /*G4cout<<"New weight brem "<<new_weight<<std::endl;
258  G4cout<<"Weight correction brem "<<w_corr<<std::endl;*/
259  fParticleChange->SetParentWeightByProcess(false);
260  fParticleChange->SetSecondaryWeightByProcess(false);
261  fParticleChange->ProposeParentWeight(new_weight);
262 
263  //Kinematic
264  //---------
266  G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
267  G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
268  G4double projectileP = std::sqrt(projectileP2);
269 
270 
271  //Use the angular model of the forward model to generate the gamma direction
272  //---------------------------------------------------------------------------
273 //Dum dynamic particle to use the model
274  G4DynamicParticle * aDynPart = new G4DynamicParticle(G4Electron::Electron(),G4ThreeVector(0.,0.,1.)*projectileP);
275 
276  //Get the element from the direct model
278  projectileKinEnergy,currentTcutForDirectSecond);
279  G4int Z=elm->GetZasInt();
280  G4double energy = aDynPart->GetTotalEnergy()-gammaEnergy;
281  G4ThreeVector projectileMomentum =
282  theDirectEMModel->GetAngularDistribution()->SampleDirection(aDynPart,energy,Z,currentMaterial)*projectileP;
283  G4double phi = projectileMomentum.getPhi();
284 
285 /*
286  //Angle of the gamma direction with the projectile taken from G4eBremsstrahlungModel
287  //------------------------------------------------
288  G4double u;
289  const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
290 
291  if (9./(9.+d) > G4UniformRand()) u = - std::log(G4UniformRand()*G4UniformRand())/a1;
292  else u = - std::log(G4UniformRand()*G4UniformRand())/a2;
293 
294  G4double theta = u*electron_mass_c2/projectileTotalEnergy;
295 
296  G4double sint = std::sin(theta);
297  G4double cost = std::cos(theta);
298 
299  G4double phi = twopi * G4UniformRand() ;
300  G4ThreeVector projectileMomentum;
301  projectileMomentum=G4ThreeVector(std::cos(phi)*sint,std::sin(phi)*sint,cost)*projectileP; //gamma frame
302 */
303  if (IsScatProjToProjCase) {//the adjoint primary is the scattered e-
304  G4ThreeVector gammaMomentum = (projectileTotalEnergy-adjointPrimTotalEnergy)*G4ThreeVector(0.,0.,1.);
305  G4ThreeVector dirProd=projectileMomentum-gammaMomentum;
306  G4double cost1 = std::cos(dirProd.angle(projectileMomentum));
307  G4double sint1 = std::sqrt(1.-cost1*cost1);
308  projectileMomentum=G4ThreeVector(std::cos(phi)*sint1,std::sin(phi)*sint1,cost1)*projectileP;
309  }
310 
311  projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection());
312 
313  if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
314  fParticleChange->ProposeTrackStatus(fStopAndKill);
315  fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
316  }
317  else {
318  fParticleChange->ProposeEnergy(projectileKinEnergy);
319  fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
320  }
321 }
323 //
325  G4double kinEnergyProj, // kinetic energy of the primary particle before the interaction
326  G4double kinEnergyProd // kinetic energy of the secondary particle
327  )
331  }
332 /*
333  return DiffCrossSectionPerVolumePrimToSecondApproximated2(aMaterial,
334  kinEnergyProj,
335  kinEnergyProd);
336 */
338  kinEnergyProj,
339  kinEnergyProd);
340 }
341 
343 //
345  const G4Material* aMaterial,
346  G4double kinEnergyProj, // kinetic energy of the primary particle before the interaction
347  G4double kinEnergyProd // kinetic energy of the secondary particle
348  )
349 {
350  G4double dCrossEprod=0.;
351  G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
352  G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
353 
354 
355  //In this approximation we consider that the secondary gammas are sampled with 1/Egamma energy distribution
356  //This is what is applied in the discrete standard model before the rejection test that make a correction
357  //The application of the same rejection function is not possible here.
358  //The differentiation of the CS over Ecut does not produce neither a good differential CS. That is due to the
359  // fact that in the discrete model the differential CS and the integrated CS are both fitted but separatly and
360  // therefore do not allow a correct numerical differentiation of the integrated CS to get the differential one.
361  // In the future we plan to use the brem secondary spectra from the G4Penelope implementation
362 
363  if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){
365  dCrossEprod=sigma/kinEnergyProd/std::log(kinEnergyProj/keV);
366  }
367  return dCrossEprod;
368 
369 }
370 
372 //
374  const G4Material* material,
375  G4double kinEnergyProj, // kinetic energy of the primary particle before the interaction
376  G4double kinEnergyProd // kinetic energy of the secondary particle
377  )
378 {
379  //In this approximation we derive the direct cross section over Tcut=gamma energy, en after apply the Migdla correction factor
380  //used in the direct model
381 
382  G4double dCrossEprod=0.;
383 
384  const G4ElementVector* theElementVector = material->GetElementVector();
385  const double* theAtomNumDensityVector = material->GetAtomicNumDensityVector();
386  G4double dum=0.;
387  G4double E1=kinEnergyProd,E2=kinEnergyProd*1.001;
388  G4double dE=E2-E1;
389  for (size_t i=0; i<material->GetNumberOfElements(); i++) {
390  G4double C1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,(*theElementVector)[i]->GetZ(),dum ,E1);
391  G4double C2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,(*theElementVector)[i]->GetZ(),dum,E2);
392  dCrossEprod += theAtomNumDensityVector[i] * (C1-C2)/dE;
393  }
394 
395  return dCrossEprod;
396 
397 }
399 //
401  G4double primEnergy,
402  G4bool IsScatProjToProjCase)
403 { if (!isDirectModelInitialised) {
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= 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 
const double C2
static G4AdjointGamma * AdjointGamma()
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:257
virtual G4double DiffCrossSectionPerVolumePrimToSecond(const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyProd)
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase(G4double PrimAdjEnergy)
const double C1
std::vector< G4Element * > G4ElementVector
virtual void SampleSecondaries(const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
G4double GetKineticEnergy() const
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
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:616
void AddEmModel(G4int, G4VEmModel *, G4VEmFluctuationModel *, const G4Region *)
static G4AdjointElectron * AdjointElectron()
G4Material * currentMaterial
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
int G4int
Definition: G4Types.hh:78
virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
static constexpr double twopi
Definition: G4SIunits.hh:76
virtual G4double GetAdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
void SetUseMatrixPerElement(G4bool aBool)
void SetSecondaryWeightByProcess(G4bool)
void SetParentWeightByProcess(G4bool)
static const G4double dE
virtual G4double DiffCrossSectionPerVolumePrimToSecond(const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyProd)
G4VEmModel * theDirectEMModel
#define G4UniformRand()
Definition: Randomize.hh:97
virtual G4ThreeVector & SampleDirection(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, const G4Material *)=0
virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
void RapidSampleSecondaries(const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)
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:321
const G4double * GetAtomicNumDensityVector() const
Definition: G4Material.hh:216
G4MaterialCutsCouple * currentCouple
void DefineCurrentMaterial(const G4MaterialCutsCouple *couple)
G4bool correct_weight_for_post_step_in_model
const G4DataVector * Initialise(const G4ParticleDefinition *part, const G4ParticleDefinition *secPart, G4double minSubRange, G4int verb)
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
G4double DiffCrossSectionPerVolumePrimToSecondApproximated2(const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyProd)
G4double GetPDGMass() const
G4double additional_weight_correction_factor_for_post_step_outside_model
G4double energy(const ThreeVector &p, const G4double m)
virtual G4double GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy, G4double Tcut=0)
G4int GetZasInt() const
Definition: G4Element.hh:132
void ProposeEnergy(G4double finalEnergy)
static constexpr double GeV
Definition: G4SIunits.hh:217
void AddSecondary(G4Track *aSecondary)
G4double GetWeight() const
static const G4double Emin
G4double DiffCrossSectionPerVolumePrimToSecondApproximated1(const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyProd)
static G4Electron * Electron()
Definition: G4Electron.cc:94
static const G4double Emax
static constexpr double MeV
Definition: G4SIunits.hh:214
G4ParticleDefinition * theAdjEquivOfDirectSecondPartDef
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
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 constexpr double keV
Definition: G4SIunits.hh:216
static G4AdjointCSManager * GetAdjointCSManager()
virtual void CorrectPostStepWeight(G4ParticleChange *fParticleChange, G4double old_weight, G4double adjointPrimKinEnergy, G4double projectileKinEnergy, G4bool IsScatProjToProjCase)
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:541
const G4Material * GetMaterial() const