Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4AdjointhIonisationModel Class Reference

#include <G4AdjointhIonisationModel.hh>

Inheritance diagram for G4AdjointhIonisationModel:
Collaboration diagram for G4AdjointhIonisationModel:

Public Member Functions

 G4AdjointhIonisationModel (G4ParticleDefinition *projectileDefinition)
 
virtual ~G4AdjointhIonisationModel ()
 
virtual void SampleSecondaries (const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)
 
void RapidSampleSecondaries (const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)
 
virtual G4double DiffCrossSectionPerAtomPrimToSecond (G4double kinEnergyProj, G4double kinEnergyProd, G4double Z, G4double A=0.)
 
virtual G4double AdjointCrossSection (const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
 
virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase (G4double PrimAdjEnergy)
 
virtual G4double GetSecondAdjEnergyMinForScatProjToProjCase (G4double PrimAdjEnergy, G4double Tcut=0)
 
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase (G4double PrimAdjEnergy)
 
virtual G4double GetSecondAdjEnergyMinForProdToProjCase (G4double PrimAdjEnergy)
 
- Public Member Functions inherited from G4VEmAdjointModel
 G4VEmAdjointModel (const G4String &nam)
 
virtual ~G4VEmAdjointModel ()
 
virtual G4double GetAdjointCrossSection (const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
 
virtual G4double DiffCrossSectionPerAtomPrimToScatPrim (G4double kinEnergyProj, G4double kinEnergyScatProj, G4double Z, G4double A=0.)
 
virtual G4double DiffCrossSectionPerVolumePrimToSecond (const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyProd)
 
virtual G4double DiffCrossSectionPerVolumePrimToScatPrim (const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyScatProj)
 
void DefineCurrentMaterial (const G4MaterialCutsCouple *couple)
 
std::vector< std::vector
< double > * > 
ComputeAdjointCrossSectionVectorPerAtomForSecond (G4double kinEnergyProd, G4double Z, G4double A=0., G4int nbin_pro_decade=10)
 
std::vector< std::vector
< double > * > 
ComputeAdjointCrossSectionVectorPerAtomForScatProj (G4double kinEnergyProd, G4double Z, G4double A=0., G4int nbin_pro_decade=10)
 
std::vector< std::vector
< double > * > 
ComputeAdjointCrossSectionVectorPerVolumeForSecond (G4Material *aMaterial, G4double kinEnergyProd, G4int nbin_pro_decade=10)
 
std::vector< std::vector
< double > * > 
ComputeAdjointCrossSectionVectorPerVolumeForScatProj (G4Material *aMaterial, G4double kinEnergyProd, G4int nbin_pro_decade=10)
 
void SetCSMatrices (std::vector< G4AdjointCSMatrix * > *Vec1CSMatrix, std::vector< G4AdjointCSMatrix * > *Vec2CSMatrix)
 
G4ParticleDefinitionGetAdjointEquivalentOfDirectPrimaryParticleDefinition ()
 
G4ParticleDefinitionGetAdjointEquivalentOfDirectSecondaryParticleDefinition ()
 
G4double GetHighEnergyLimit ()
 
G4double GetLowEnergyLimit ()
 
void SetHighEnergyLimit (G4double aVal)
 
void SetLowEnergyLimit (G4double aVal)
 
void DefineDirectEMModel (G4VEmModel *aModel)
 
void SetAdjointEquivalentOfDirectPrimaryParticleDefinition (G4ParticleDefinition *aPart)
 
void SetAdjointEquivalentOfDirectSecondaryParticleDefinition (G4ParticleDefinition *aPart)
 
void SetSecondPartOfSameType (G4bool aBool)
 
G4bool GetSecondPartOfSameType ()
 
void SetUseMatrix (G4bool aBool)
 
void SetUseMatrixPerElement (G4bool aBool)
 
void SetUseOnlyOneMatrixForAllElements (G4bool aBool)
 
void SetApplyCutInRange (G4bool aBool)
 
G4bool GetUseMatrix ()
 
G4bool GetUseMatrixPerElement ()
 
G4bool GetUseOnlyOneMatrixForAllElements ()
 
G4bool GetApplyCutInRange ()
 
G4String GetName ()
 
virtual void SetCSBiasingFactor (G4double aVal)
 
void SetCorrectWeightForPostStepInModel (G4bool aBool)
 
void SetAdditionalWeightCorrectionFactorForPostStepOutsideModel (G4double factor)
 

Additional Inherited Members

- Protected Member Functions inherited from G4VEmAdjointModel
G4double DiffCrossSectionFunction1 (G4double kinEnergyProj)
 
G4double DiffCrossSectionFunction2 (G4double kinEnergyProj)
 
G4double DiffCrossSectionPerVolumeFunctionForIntegrationOverEkinProj (G4double EkinProd)
 
G4double SampleAdjSecEnergyFromCSMatrix (size_t MatrixIndex, G4double prim_energy, G4bool IsScatProjToProjCase)
 
G4double SampleAdjSecEnergyFromCSMatrix (G4double prim_energy, G4bool IsScatProjToProjCase)
 
void SelectCSMatrix (G4bool IsScatProjToProjCase)
 
virtual G4double SampleAdjSecEnergyFromDiffCrossSectionPerAtom (G4double prim_energy, G4bool IsScatProjToProjCase)
 
virtual void CorrectPostStepWeight (G4ParticleChange *fParticleChange, G4double old_weight, G4double adjointPrimKinEnergy, G4double projectileKinEnergy, G4bool IsScatProjToProjCase)
 
- Protected Attributes inherited from G4VEmAdjointModel
G4VEmModeltheDirectEMModel
 
G4VParticleChangepParticleChange
 
const G4String name
 
G4int ASelectedNucleus
 
G4int ZSelectedNucleus
 
G4MaterialSelectedMaterial
 
G4double kinEnergyProdForIntegration
 
G4double kinEnergyScatProjForIntegration
 
G4double kinEnergyProjForIntegration
 
std::vector< G4AdjointCSMatrix * > * pOnCSMatrixForProdToProjBackwardScattering
 
std::vector< G4AdjointCSMatrix * > * pOnCSMatrixForScatProjToProjBackwardScattering
 
std::vector< G4doubleCS_Vs_ElementForScatProjToProjCase
 
std::vector< G4doubleCS_Vs_ElementForProdToProjCase
 
G4double lastCS
 
G4double lastAdjointCSForScatProjToProjCase
 
G4double lastAdjointCSForProdToProjCase
 
G4ParticleDefinitiontheAdjEquivOfDirectPrimPartDef
 
G4ParticleDefinitiontheAdjEquivOfDirectSecondPartDef
 
G4ParticleDefinitiontheDirectPrimaryPartDef
 
G4bool second_part_of_same_type
 
G4double preStepEnergy
 
G4MaterialcurrentMaterial
 
G4MaterialCutsCouplecurrentCouple
 
size_t currentMaterialIndex
 
size_t currentCoupleIndex
 
G4double currentTcutForDirectPrim
 
G4double currentTcutForDirectSecond
 
G4bool ApplyCutInRange
 
G4double mass_ratio_product
 
G4double mass_ratio_projectile
 
G4double HighEnergyLimit
 
G4double LowEnergyLimit
 
G4double CS_biasing_factor
 
G4bool UseMatrix
 
G4bool UseMatrixPerElement
 
G4bool UseOnlyOneMatrixForAllElements
 
size_t indexOfUsedCrossSectionMatrix
 
size_t model_index
 
G4bool correct_weight_for_post_step_in_model
 
G4double additional_weight_correction_factor_for_post_step_outside_model
 

Detailed Description

Definition at line 72 of file G4AdjointhIonisationModel.hh.

Constructor & Destructor Documentation

G4AdjointhIonisationModel::G4AdjointhIonisationModel ( G4ParticleDefinition projectileDefinition)

Definition at line 46 of file G4AdjointhIonisationModel.cc.

46  :
47  G4VEmAdjointModel("Adjoint_hIonisation")
48 {
49 
50 
51 
52  UseMatrix =true;
53  UseMatrixPerElement = true;
54  ApplyCutInRange = true;
58 
59  //The direct EM Modfel is taken has BetheBloch it is only used for the computation
60  // of the differential cross section.
61  //The Bragg model could be used as an alternative as it offers the same differential cross section
62 
63  theDirectEMModel = new G4BetheBlochModel(projectileDefinition);
64  theBraggDirectEMModel = new G4BraggModel(projectileDefinition);
66 
67  theDirectPrimaryPartDef = projectileDefinition;
69  if (projectileDefinition == G4Proton::Proton()) {
71 
72  }
73 
74  DefineProjectileProperty();
75 }
G4ParticleDefinition * theDirectPrimaryPartDef
static G4AdjointElectron * AdjointElectron()
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
G4VEmModel * theDirectEMModel
G4VEmAdjointModel(const G4String &nam)
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4bool UseOnlyOneMatrixForAllElements
static G4AdjointProton * AdjointProton()
G4ParticleDefinition * theAdjEquivOfDirectSecondPartDef

Here is the call graph for this function:

G4AdjointhIonisationModel::~G4AdjointhIonisationModel ( )
virtual

Definition at line 78 of file G4AdjointhIonisationModel.cc.

79 {;}

Member Function Documentation

G4double G4AdjointhIonisationModel::AdjointCrossSection ( const G4MaterialCutsCouple aCouple,
G4double  primEnergy,
G4bool  IsScatProjToProjCase 
)
virtual

Reimplemented from G4VEmAdjointModel.

Definition at line 438 of file G4AdjointhIonisationModel.cc.

441 {
442  if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase);
443  DefineCurrentMaterial(aCouple);
444 
445 
447 
448  if (!IsScatProjToProjCase ){
449  G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy);
450  G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy);
451  if (Emax_proj>Emin_proj && primEnergy > currentTcutForDirectSecond) {
452  Cross*=(1./Emin_proj -1./Emax_proj)/primEnergy;
453  }
454  else Cross=0.;
455 
456 
457 
458 
459 
460 
461  }
462  else {
465  G4double diff1=Emin_proj-primEnergy;
466  G4double diff2=Emax_proj-primEnergy;
467  G4double t1=(1./diff1+1./Emin_proj-1./diff2-1./Emax_proj)/primEnergy;
468  //G4double t2=2.*std::log(diff2*Emin_proj/Emax_proj/diff1)/primEnergy/primEnergy;
469  G4double t2=2.*std::log(Emax_proj/Emin_proj)/primEnergy/primEnergy;
470  Cross*=(t1+t2);
471 
472 
473  }
474  lastCS =Cross;
475  return Cross;
476 }
virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
G4Material * currentMaterial
virtual G4double GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy, G4double Tcut=0)
G4double GetElectronDensity() const
Definition: G4Material.hh:217
void DefineCurrentMaterial(const G4MaterialCutsCouple *couple)
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
tuple t1
Definition: plottest35.py:33
G4double currentTcutForDirectSecond
double G4double
Definition: G4Types.hh:76
virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase(G4double PrimAdjEnergy)

Here is the call graph for this function:

G4double G4AdjointhIonisationModel::DiffCrossSectionPerAtomPrimToSecond ( G4double  kinEnergyProj,
G4double  kinEnergyProd,
G4double  Z,
G4double  A = 0. 
)
virtual

Reimplemented from G4VEmAdjointModel.

Definition at line 306 of file G4AdjointhIonisationModel.cc.

311 {//Probably that here the Bragg Model should be also used for kinEnergyProj/nuc<2MeV
312 
313 
314 
315  G4double dSigmadEprod=0;
316  G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
317  G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
318 
319 
320  if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ //the produced particle should have a kinetic energy smaller than the projectile
321  G4double Tmax=kinEnergyProj;
322 
323  G4double E1=kinEnergyProd;
324  G4double E2=kinEnergyProd*1.000001;
325  G4double dE=(E2-E1);
326  G4double sigma1,sigma2;
327  if (kinEnergyProj >2.*MeV){
328  sigma1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
329  sigma2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
330  }
331  else {
332  sigma1=theBraggDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
333  sigma2=theBraggDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
334  }
335 
336 
337  dSigmadEprod=(sigma1-sigma2)/dE;
338  if (dSigmadEprod>1.) {
339  G4cout<<"sigma1 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma1<<G4endl;
340  G4cout<<"sigma2 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma2<<G4endl;
341  G4cout<<"dsigma "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<dSigmadEprod<<G4endl;
342 
343  }
344 
345 
346 
347  //correction of differential cross section at high energy to correct for the suppression of particle at secondary at high
348  //energy used in the Bethe Bloch Model. This correction consist to multiply by g the probability function used
349  //to test the rejection of a secondary
350  //-------------------------
351 
352  //Source code taken from G4BetheBlochModel::SampleSecondaries
353 
354  G4double deltaKinEnergy = kinEnergyProd;
355 
356  //Part of the taken code
357  //----------------------
358 
359 
360 
361  // projectile formfactor - suppresion of high energy
362  // delta-electron production at high energy
363  G4double x = formfact*deltaKinEnergy;
364  if(x > 1.e-6) {
365 
366 
367  G4double totEnergy = kinEnergyProj + mass;
368  G4double etot2 = totEnergy*totEnergy;
369  G4double beta2 = kinEnergyProj*(kinEnergyProj + 2.0*mass)/etot2;
370  G4double f;
371  G4double f1 = 0.0;
372  f = 1.0 - beta2*deltaKinEnergy/Tmax;
373  if( 0.5 == spin ) {
374  f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
375  f += f1;
376  }
377  G4double x1 = 1.0 + x;
378  G4double gg = 1.0/(x1*x1);
379  if( 0.5 == spin ) {
380  G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
381  gg *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
382  }
383  if(gg > 1.0) {
384  G4cout << "### G4BetheBlochModel in Adjoint Sim WARNING: g= " << g
385  << G4endl;
386  gg=1.;
387  }
388  //G4cout<<"gg"<<gg<<G4endl;
389  dSigmadEprod*=gg;
390  }
391 
392  }
393 
394  return dSigmadEprod;
395 }
G4ParticleDefinition * theDirectPrimaryPartDef
tuple x
Definition: test.py:50
function g(Y1, Y2, PT2)
Definition: hijing1.383.f:5205
static const G4double dE
G4VEmModel * theDirectEMModel
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
float electron_mass_c2
Definition: hepunit.py:274
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:321
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
double G4double
Definition: G4Types.hh:76
virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase(G4double PrimAdjEnergy)

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMaxForProdToProjCase ( G4double  PrimAdjEnergy)
virtual

Reimplemented from G4VEmAdjointModel.

Definition at line 491 of file G4AdjointhIonisationModel.cc.

492 { return HighEnergyLimit;
493 }

Here is the caller graph for this function:

G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMaxForScatProjToProjCase ( G4double  PrimAdjEnergy)
virtual

Reimplemented from G4VEmAdjointModel.

Definition at line 479 of file G4AdjointhIonisationModel.cc.

480 {
481  G4double Tmax=PrimAdjEnergy*one_plus_ratio_2/(one_minus_ratio_2-2.*ratio*PrimAdjEnergy/mass);
482  return Tmax;
483 }
double G4double
Definition: G4Types.hh:76

Here is the caller graph for this function:

G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMinForProdToProjCase ( G4double  PrimAdjEnergy)
virtual

Reimplemented from G4VEmAdjointModel.

Definition at line 496 of file G4AdjointhIonisationModel.cc.

497 { G4double Tmin= (2*PrimAdjEnergy-4*mass + std::sqrt(4.*PrimAdjEnergy*PrimAdjEnergy +16.*mass*mass + 8.*PrimAdjEnergy*mass*(1/ratio +ratio)))/4.;
498  return Tmin;
499 }
double G4double
Definition: G4Types.hh:76

Here is the caller graph for this function:

G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMinForScatProjToProjCase ( G4double  PrimAdjEnergy,
G4double  Tcut = 0 
)
virtual

Reimplemented from G4VEmAdjointModel.

Definition at line 486 of file G4AdjointhIonisationModel.cc.

487 { return PrimAdjEnergy+Tcut;
488 }

Here is the caller graph for this function:

void G4AdjointhIonisationModel::RapidSampleSecondaries ( const G4Track aTrack,
G4bool  IsScatProjToProjCase,
G4ParticleChange fParticleChange 
)

Definition at line 160 of file G4AdjointhIonisationModel.cc.

163 {
164 
165  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
167 
168 
169  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
170  G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
171 
172  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
173  return;
174  }
175 
176  G4double projectileKinEnergy =0.;
177  G4double eEnergy=0.;
179  if (!IsScatProjToProjCase){//1/E^2 distribution
180 
181  eEnergy=adjointPrimKinEnergy;
182  G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
184  if (Emin>=Emax) return;
185  G4double a=1./Emax;
186  G4double b=1./Emin;
187  newCS=newCS*(b-a)/eEnergy;
188 
189  projectileKinEnergy =1./(b- (b-a)*G4UniformRand());
190 
191 
192  }
193  else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
195  if (Emin>=Emax) return;
196  G4double diff1=Emin-adjointPrimKinEnergy;
197  G4double diff2=Emax-adjointPrimKinEnergy;
198 
199  G4double t1=adjointPrimKinEnergy*(1./diff1-1./diff2);
200  G4double t2=adjointPrimKinEnergy*(1./Emin-1./Emax);
201  /*G4double f31=diff1/Emin;
202  G4double f32=diff2/Emax/f31;*/
203  G4double t3=2.*std::log(Emax/Emin);
204  G4double sum_t=t1+t2+t3;
205  newCS=newCS*sum_t/adjointPrimKinEnergy/adjointPrimKinEnergy;
206  G4double t=G4UniformRand()*sum_t;
207  if (t <=t1 ){
208  G4double q= G4UniformRand()*t1/adjointPrimKinEnergy ;
209  projectileKinEnergy =adjointPrimKinEnergy +1./(1./diff1-q);
210 
211  }
212  else if (t <=t2 ) {
213  G4double q= G4UniformRand()*t2/adjointPrimKinEnergy;
214  projectileKinEnergy =1./(1./Emin-q);
215  }
216  else {
217  projectileKinEnergy=Emin*std::pow(Emax/Emin,G4UniformRand());
218 
219  }
220  eEnergy=projectileKinEnergy-adjointPrimKinEnergy;
221 
222 
223  }
224 
225 
226 
227  G4double diffCS_perAtom_Used=twopi_mc2_rcl2*mass*adjointPrimKinEnergy/projectileKinEnergy/projectileKinEnergy/eEnergy/eEnergy;
228 
229 
230 
231  //Weight correction
232  //-----------------------
233  //First w_corr is set to the ratio between adjoint total CS and fwd total CS
235 
236  //G4cout<<w_corr<<G4endl;
237  w_corr*=newCS/lastCS;
238  //G4cout<<w_corr<<G4endl;
239  //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
240  //Here we consider the true diffCS as the one obtained by the numerical differentiation over Tcut of the direct CS
241 
242  G4double diffCS = DiffCrossSectionPerAtomPrimToSecond(projectileKinEnergy, eEnergy,1,1);
243  w_corr*=diffCS/diffCS_perAtom_Used;
244  //G4cout<<w_corr<<G4endl;
245 
246  G4double new_weight = aTrack.GetWeight()*w_corr;
247  fParticleChange->SetParentWeightByProcess(false);
248  fParticleChange->SetSecondaryWeightByProcess(false);
249  fParticleChange->ProposeParentWeight(new_weight);
250 
251 
252 
253 
254  //Kinematic:
255  //we consider a two body elastic scattering for the forward processes where the projectile knock on an e- at rest and gives
256  // him part of its energy
257  //----------------------------------------------------------------------------------------
258 
260  G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
261  G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
262 
263 
264 
265  //Companion
266  //-----------
268  if (IsScatProjToProjCase) {
270  }
271  G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
272  G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;
273 
274 
275  //Projectile momentum
276  //--------------------
277  G4double P_parallel = (adjointPrimP*adjointPrimP + projectileP2 - companionP2)/(2.*adjointPrimP);
278  G4double P_perp = std::sqrt( projectileP2 - P_parallel*P_parallel);
279  G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
280  G4double phi =G4UniformRand()*2.*3.1415926;
281  G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
282  projectileMomentum.rotateUz(dir_parallel);
283 
284 
285 
286  if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
287  fParticleChange->ProposeTrackStatus(fStopAndKill);
288  fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
289  //G4cout<<"projectileMomentum "<<projectileMomentum<<G4endl;
290  }
291  else {
292  fParticleChange->ProposeEnergy(projectileKinEnergy);
293  fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
294  }
295 
296 
297 
298 
299 
300 
301 
302 }
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4double GetPostStepWeightCorrection()
const G4DynamicParticle * GetDynamicParticle() const
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
void ProposeParentWeight(G4double finalWeight)
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
G4Material * currentMaterial
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
virtual G4double GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy, G4double Tcut=0)
G4double GetTotalMomentum() const
tuple b
Definition: test.py:12
void SetSecondaryWeightByProcess(G4bool)
void SetParentWeightByProcess(G4bool)
#define G4UniformRand()
Definition: Randomize.hh:97
G4double GetElectronDensity() const
Definition: G4Material.hh:217
const G4ThreeVector & GetMomentumDirection() const
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
void DefineCurrentMaterial(const G4MaterialCutsCouple *couple)
G4double GetPDGMass() const
Hep3Vector unit() const
tuple t1
Definition: plottest35.py:33
void ProposeEnergy(G4double finalEnergy)
void AddSecondary(G4Track *aSecondary)
G4double GetWeight() const
static const G4double Emin
static const G4double Emax
G4ParticleDefinition * theAdjEquivOfDirectSecondPartDef
G4double currentTcutForDirectSecond
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
static G4AdjointCSManager * GetAdjointCSManager()
virtual G4double DiffCrossSectionPerAtomPrimToSecond(G4double kinEnergyProj, G4double kinEnergyProd, G4double Z, G4double A=0.)
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase(G4double PrimAdjEnergy)

Here is the call graph for this function:

Here is the caller graph for this function:

void G4AdjointhIonisationModel::SampleSecondaries ( const G4Track aTrack,
G4bool  IsScatProjToProjCase,
G4ParticleChange fParticleChange 
)
virtual

Implements G4VEmAdjointModel.

Definition at line 84 of file G4AdjointhIonisationModel.cc.

87 {
88  if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange);
89 
90  const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
91 
92  //Elastic inverse scattering
93  //---------------------------------------------------------
94  G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
95  G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
96 
97  if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
98  return;
99  }
100 
101  //Sample secondary energy
102  //-----------------------
103  G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, IsScatProjToProjCase);
104  CorrectPostStepWeight(fParticleChange,
105  aTrack.GetWeight(),
106  adjointPrimKinEnergy,
107  projectileKinEnergy,
108  IsScatProjToProjCase); //Caution !!!this weight correction should be always applied
109 
110 
111  //Kinematic:
112  //we consider a two body elastic scattering for the forward processes where the projectile knock on an e- at rest and gives
113  // him part of its energy
114  //----------------------------------------------------------------------------------------
115 
117  G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
118  G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
119 
120 
121 
122  //Companion
123  //-----------
125  if (IsScatProjToProjCase) {
127  }
128  G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
129  G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;
130 
131 
132  //Projectile momentum
133  //--------------------
134  G4double P_parallel = (adjointPrimP*adjointPrimP + projectileP2 - companionP2)/(2.*adjointPrimP);
135  G4double P_perp = std::sqrt( projectileP2 - P_parallel*P_parallel);
136  G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
137  G4double phi =G4UniformRand()*2.*3.1415926;
138  G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
139  projectileMomentum.rotateUz(dir_parallel);
140 
141 
142 
143  if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
144  fParticleChange->ProposeTrackStatus(fStopAndKill);
145  fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
146  //G4cout<<"projectileMomentum "<<projectileMomentum<<G4endl;
147  }
148  else {
149  fParticleChange->ProposeEnergy(projectileKinEnergy);
150  fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
151  }
152 
153 
154 
155 
156 }
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
const G4DynamicParticle * GetDynamicParticle() const
G4double SampleAdjSecEnergyFromCSMatrix(size_t MatrixIndex, G4double prim_energy, G4bool IsScatProjToProjCase)
void RapidSampleSecondaries(const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
G4double GetTotalMomentum() const
#define G4UniformRand()
Definition: Randomize.hh:97
const G4ThreeVector & GetMomentumDirection() const
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
G4double GetPDGMass() const
Hep3Vector unit() const
void ProposeEnergy(G4double finalEnergy)
void AddSecondary(G4Track *aSecondary)
G4double GetWeight() const
G4ParticleDefinition * theAdjEquivOfDirectSecondPartDef
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
virtual void CorrectPostStepWeight(G4ParticleChange *fParticleChange, G4double old_weight, G4double adjointPrimKinEnergy, G4double projectileKinEnergy, G4bool IsScatProjToProjCase)

Here is the call graph for this function:


The documentation for this class was generated from the following files: