lars/geant4
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Import Geant4 9.3.0 source tree

Browse Source
This commit is contained in:
Gabriele Cosmo
2016-06-09 16:15:05 +02:00
parent b79225fb37
commit 74cad5e589
3877 changed files with 234205 additions and 167127 deletions
Download Patch File Download Diff File Expand all files Collapse all files
source/processes/electromagnetic/adjoint/src/G4AdjointAlongStepWeightCorrection.cc
+33 -11
View File
@@ -23,13 +23,15 @@
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
// $Id: G4AdjointAlongStepWeightCorrection.cc,v 1.5 2009/11/23 09:02:35 gcosmo Exp $
// GEANT4 tag $Name: geant4-09-03 $
//
#include "G4AdjointAlongStepWeightCorrection.hh"
#include "G4Step.hh"
#include "G4ParticleDefinition.hh"
#include "G4VParticleChange.hh"
#include "G4AdjointCSManager.hh"
///////////////////////////////////////////////////////
//
@@ -44,28 +46,19 @@ G4AdjointAlongStepWeightCorrection::G4AdjointAlongStepWeightCorrection(const G4S
G4AdjointAlongStepWeightCorrection::~G4AdjointAlongStepWeightCorrection()
{;
}
///////////////////////////////////////////////////////
//
void G4AdjointAlongStepWeightCorrection::PreparePhysicsTable(
const G4ParticleDefinition& )
{
;
}
///////////////////////////////////////////////////////
//
void G4AdjointAlongStepWeightCorrection::BuildPhysicsTable(const G4ParticleDefinition& )
{;
}
///////////////////////////////////////////////////////
//
G4VParticleChange* G4AdjointAlongStepWeightCorrection::AlongStepDoIt(const G4Track& track,
@@ -84,9 +77,28 @@ G4VParticleChange* G4AdjointAlongStepWeightCorrection::AlongStepDoIt(const G4Tra
G4double weight_correction=G4AdjointCSManager::GetAdjointCSManager()->GetContinuousWeightCorrection(thePartDef,
preStepKinEnergy,Tkin, currentCouple,length);
G4double new_weight=weight_correction*track.GetWeight();
//if (weight_correction >2.) new_weight=1.e-300;
//The following test check for zero weight.
//This happens after weight correction of gamma for photo electric effect.
//When the new weight is 0 it will be later on consider as nan by G4.
//Therefore we do put a lower limit of 1.e-300. for new_weight
//Correction by L.Desorgher on 15 July 2009
#ifdef WIN32
if (!!_isnan(new_weight) || new_weight==0){
#else
if (std::isnan(new_weight) || new_weight==0){
#endif
//G4cout<<new_weight<<'\t'<<weight_correction<<'\t'<<track.GetWeight()<<G4endl;
new_weight=1.e-300;
}
//G4cout<<new_weight<<'\t'<<weight_correction<<'\t'<<track.GetWeight()<<G4endl;
fParticleChange->SetParentWeightByProcess(false);
fParticleChange->SetSecondaryWeightByProcess(false);
fParticleChange->ProposeParentWeight(new_weight);
@@ -95,3 +107,13 @@ G4VParticleChange* G4AdjointAlongStepWeightCorrection::AlongStepDoIt(const G4Tra
return fParticleChange;
}
///////////////////////////////////////////////////////
//
G4double G4AdjointAlongStepWeightCorrection::GetContinuousStepLimit(const G4Track& track,
G4double , G4double , G4double& )
{
G4double x = DBL_MAX;
DefineMaterial(track.GetMaterialCutsCouple());
preStepKinEnergy = track.GetKineticEnergy();
return x;
}
source/processes/electromagnetic/adjoint/src/G4AdjointBremsstrahlungModel.cc
+274 -544
View File
@@ -23,467 +23,59 @@
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
// $Id: G4AdjointBremsstrahlungModel.cc,v 1.5 2009/12/16 17:50:01 gunter Exp $
// GEANT4 tag $Name: geant4-09-03 $
//
#include "G4AdjointBremsstrahlungModel.hh"
#include "G4AdjointCSManager.hh"
#include "G4Integrator.hh"
#include "G4TrackStatus.hh"
#include "G4ParticleChange.hh"
#include "G4AdjointElectron.hh"
#include "G4AdjointGamma.hh"
#include "G4Electron.hh"
#include "G4Timer.hh"
//#include "G4PenelopeBremsstrahlungModel.hh"
////////////////////////////////////////////////////////////////////////////////
//
G4AdjointBremsstrahlungModel::G4AdjointBremsstrahlungModel():
G4VEmAdjointModel("AdjointBremModel"),
probsup(1.0),
MigdalConstant(classic_electr_radius*electron_Compton_length*electron_Compton_length/pi),
LPMconstant(fine_structure_const*electron_mass_c2*electron_mass_c2/(4.*pi*hbarc)),
theLPMflag(true)
{ isElectron= true;
SetUseMatrix(true);
G4VEmAdjointModel("AdjointeBremModel"),
MigdalConstant(classic_electr_radius*electron_Compton_length*electron_Compton_length*4.0*pi)
{
SetUseMatrix(false);
SetUseMatrixPerElement(false);
theDirectStdBremModel = new G4eBremsstrahlungModel(G4Electron::Electron(),"TheDirecteBremModel");
theDirectEMModel=theDirectStdBremModel;
// theDirectPenelopeBremModel =0;
SetApplyCutInRange(true);
SetIsIonisation(false);
highKinEnergy= 100.*TeV;
lowKinEnergy = 1.0*keV;
theTimer =new G4Timer();
theTimer->Start();
InitialiseParameters();
theTimer->Stop();
G4cout<<"Time elapsed in second for the initialidation of AdjointBrem "<<theTimer->GetRealElapsed()<<std::endl;
theAdjEquivOfDirectPrimPartDef =G4AdjointElectron::AdjointElectron();
theAdjEquivOfDirectSecondPartDef=G4AdjointGamma::AdjointGamma();
theDirectPrimaryPartDef=G4Electron::Electron();
second_part_of_same_type=false;
ModeldCS="MODEL1";
/*UsePenelopeModel=false;
if (UsePenelopeModel) {
G4PenelopeBremsstrahlungModel* thePenelopeModel = new G4PenelopeBremsstrahlungModel(G4Electron::Electron(),"PenelopeBrem");
theEmModelManagerForFwdModels = new G4EmModelManager();
isPenelopeModelInitialised = false;
G4VEmFluctuationModel* f=0;
G4Region* r=0;
theDirectEMModel=thePenelopeModel;
theEmModelManagerForFwdModels->AddEmModel(1, thePenelopeModel, f, r);
}
*/
}
////////////////////////////////////////////////////////////////////////////////
//
G4AdjointBremsstrahlungModel::~G4AdjointBremsstrahlungModel()
{;}
////////////////////////////////////////////////////////////////////////////////
//
/*G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond(
const G4Material* aMaterial,
G4double kinEnergyProj, // kinetic energy of the primary particle before the interaction
G4double kinEnergyProd // kinetic energy of the secondary particle
)
{
static const G4double
ah10 = 4.67733E+00, ah11 =-6.19012E-01, ah12 = 2.02225E-02,
ah20 =-7.34101E+00, ah21 = 1.00462E+00, ah22 =-3.20985E-02,
ah30 = 2.93119E+00, ah31 =-4.03761E-01, ah32 = 1.25153E-02;
static const G4double
bh10 = 4.23071E+00, bh11 =-6.10995E-01, bh12 = 1.95531E-02,
bh20 =-7.12527E+00, bh21 = 9.69160E-01, bh22 =-2.74255E-02,
bh30 = 2.69925E+00, bh31 =-3.63283E-01, bh32 = 9.55316E-03;
static const G4double
al00 =-2.05398E+00, al01 = 2.38815E-02, al02 = 5.25483E-04,
al10 =-7.69748E-02, al11 =-6.91499E-02, al12 = 2.22453E-03,
al20 = 4.06463E-02, al21 =-1.01281E-02, al22 = 3.40919E-04;
static const G4double
bl00 = 1.04133E+00, bl01 =-9.43291E-03, bl02 =-4.54758E-04,
bl10 = 1.19253E-01, bl11 = 4.07467E-02, bl12 =-1.30718E-03,
bl20 =-1.59391E-02, bl21 = 7.27752E-03, bl22 =-1.94405E-04;
static const G4double tlow = 1.*MeV;
G4double dCrossEprod=0.;
G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){
G4double cross = 0.0;
G4double E1=kinEnergyProd;
G4double E2=kinEnergyProd*1.000000001;
G4double dE=(E2-E1);
const G4ElementVector* theElementVector = aMaterial->GetElementVector();
const G4double* theAtomNumDensityVector = aMaterial->GetAtomicNumDensityVector();
G4double dum=0.;
for (size_t i=0; i<aMaterial->GetNumberOfElements(); i++) {
G4double fac=
cross += theAtomNumDensityVector[i] * theDirectEMModel->ComputeCrossSectionPerAtom(G4Electron::Electron(),
kinEnergyProj, (*theElementVector)[i]->GetZ(), dum,E1);
}
dCrossEprod=(cross1-cross2)/dE; //first term
//Now come the correction
//-----------------------
//First compute fsig for E1
//-------------------------
G4double totalEnergy = kinEnergyProj+electron_mass_c2 ;
G4double kp2 = MigdalConstant*totalEnergy*totalEnergy
*(aMaterial->GetElectronDensity());
G4double fsig = 0.;
G4int nmax = 100;
G4double vmin=std::log(E1);
G4double vmax=std::log(kinEnergyProj) ;
G4int nn = (G4int)(nmax*(vmax-vmin)/(std::log(highKinEnergy)-vmin));
G4double u,fac,c,v,dv,y ;
if(nn > 0) {
dv = (vmax-vmin)/nn ;
v = vmin-dv ;
for(G4int n=0; n<=nn; n++) {
v += dv;
u = std::exp(v);
fac = SupressionFunction(aMaterial, kinEnergyProj, u);
y = u/kinEnergyProj;
fac *= (4.-4.*y+3.*y*y)/3.;
fac *= probsup*(u*u/(u*u+kp2))+1.-probsup;
if ((n==0)||(n==nn)) c=0.5;
else c=1. ;
fac *= c;
fsig += fac;
}
y = E1/kinEnergyProj ;
fsig *=dv/(-4.*std::log(y)/3.-4.*(1.-y)/3.+0.5*(1.-y*y));
}
else {
fsig = 1.;
}
if (fsig > 1.) fsig = 1.;
dCrossEprod*=fsig;
//return dCrossEprod;
//Now we compute dfsig
//-------------------------
G4double dfsig = 0.;
nn=20;
vmax=std::log(E2) ;
dv = (vmax-vmin)/nn ;
v = vmin-dv ;
for(G4int n=0; n<=nn; n++) {
v += dv;
u = std::exp(v);
fac = SupressionFunction(aMaterial, kinEnergyProj, u);
y = u/kinEnergyProj;
fac *= (4.-4.*y+3.*y*y)/3.;
fac *= probsup*(u*u/(u*u+kp2))+1.-probsup;
if ((n==0)||(n==nn)) c=0.5;
else c=1. ;
fac *= c;
dfsig += fac;
}
y = E1/kinEnergyProj;
dfsig *=dv/(-4.*std::log(y)/3.-4.*(1.-y)/3.+0.5*(1.-y*y));
dCrossEprod+=dfsig*cross1/dE;
}
return dCrossEprod;
}
*/
G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond(const G4Material* aMaterial,
G4double kinEnergyProj, // kinetic energy of the primary particle before the interaction
G4double kinEnergyProd // kinetic energy of the secondary particle
)
{if (ModeldCS=="MODEL2") return DiffCrossSectionPerVolumePrimToSecond2(aMaterial,
kinEnergyProj, // kinetic energy of the primary particle before the interaction
kinEnergyProd);
if (ModeldCS=="MODEL3") return DiffCrossSectionPerVolumePrimToSecond3(aMaterial,
kinEnergyProj, // kinetic energy of the primary particle before the interaction
kinEnergyProd);
return DiffCrossSectionPerVolumePrimToSecond1(aMaterial,
kinEnergyProj, // kinetic energy of the primary particle before the interaction
kinEnergyProd);
}
////////////////////////////////////////////////////////////////////////////////
// the one used till now
G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond1(
const G4Material* aMaterial,
G4double kinEnergyProj, // kinetic energy of the primary particle before the interaction
G4double kinEnergyProd // kinetic energy of the secondary particle
)
{
G4double dCrossEprod=0.;
G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){
G4double cross1 = 0.0;
G4double cross2 = 0.0;
G4double E1=kinEnergyProd;
G4double E2=kinEnergyProd*1.01;
G4double dE=(E2-E1);
const G4ElementVector* theElementVector = aMaterial->GetElementVector();
const G4double* theAtomNumDensityVector = aMaterial->GetAtomicNumDensityVector();
G4double dum=0.;
for (size_t i=0; i<aMaterial->GetNumberOfElements(); i++) {
cross1 += theAtomNumDensityVector[i] * theDirectEMModel->ComputeCrossSectionPerAtom(G4Electron::Electron(),
kinEnergyProj, (*theElementVector)[i]->GetZ(), dum,E1);
cross2 += theAtomNumDensityVector[i] * theDirectEMModel->ComputeCrossSectionPerAtom(G4Electron::Electron(),
kinEnergyProj, (*theElementVector)[i]->GetZ(), dum, E2);
}
dCrossEprod=(cross1-cross2)/dE; //first term
//Now come the correction
//-----------------------
//First compute fsig for E1
//-------------------------
G4double totalEnergy = kinEnergyProj+electron_mass_c2 ;
G4double kp2 = MigdalConstant*totalEnergy*totalEnergy
*(aMaterial->GetElectronDensity());
G4double fsig1 = 0.;
G4int nmax = 100;
G4double vmin=std::log(E1);
G4double vmax=std::log(kinEnergyProj) ;
G4int nn = (G4int)(nmax*(vmax-vmin)/(std::log(highKinEnergy)-vmin));
G4double u,fac,c,v,dv,y ;
if(nn > 0) {
dv = (vmax-vmin)/nn ;
v = vmin-dv ;
for(G4int n=0; n<=nn; n++) {
v += dv;
u = std::exp(v);
fac = SupressionFunction(aMaterial, kinEnergyProj, u);
y = u/kinEnergyProj;
fac *= (4.-4.*y+3.*y*y)/3.;
fac *= probsup*(u*u/(u*u+kp2))+1.-probsup;
if ((n==0)||(n==nn)) c=0.5;
else c=1. ;
fac *= c;
fsig1 += fac;
}
y = E1/kinEnergyProj ;
fsig1 *=dv/(-4.*std::log(y)/3.-4.*(1.-y)/3.+0.5*(1.-y*y));
}
else {
fsig1 = 1.;
}
if (fsig1 > 1.) fsig1 = 1.;
dCrossEprod*=fsig1;
G4double fsig2 = 0.;
vmin=std::log(E2);
nn = (G4int)(nmax*(vmax-vmin)/(std::log(highKinEnergy)-vmin));
if(nn > 0) {
dv = (vmax-vmin)/nn ;
v = vmin-dv ;
for(G4int n=0; n<=nn; n++) {
v += dv;
u = std::exp(v);
fac = SupressionFunction(aMaterial, kinEnergyProj, u);
y = u/kinEnergyProj;
fac *= (4.-4.*y+3.*y*y)/3.;
fac *= probsup*(u*u/(u*u+kp2))+1.-probsup;
if ((n==0)||(n==nn)) c=0.5;
else c=1. ;
fac *= c;
fsig2 += fac;
}
y = E2/kinEnergyProj ;
fsig2 *=dv/(-4.*std::log(y)/3.-4.*(1.-y)/3.+0.5*(1.-y*y));
}
else {
fsig2 = 1.;
}
if (fsig2 > 1.) fsig2 = 1.;
G4double dfsig=(fsig2-fsig1);
dCrossEprod+=dfsig*cross1/dE;
dCrossEprod=(fsig1*cross1-fsig2*cross2)/dE;
/*if (fsig < 1.){
//Now we compute dfsig
//-------------------------
G4double dfsig = 0.;
nn=20;
vmax=std::log(E2) ;
dv = (vmax-vmin)/nn ;
v = vmin-dv ;
for(G4int n=0; n<=nn; n++) {
v += dv;
u = std::exp(v);
fac = SupressionFunction(aMaterial, kinEnergyProj, u);
y = u/kinEnergyProj;
fac *= (4.-4.*y+3.*y*y)/3.;
fac *= probsup*(u*u/(u*u+kp2))+1.-probsup;
if ((n==0)||(n==nn)) c=0.5;
else c=1. ;
fac *= c;
dfsig += fac;
}
y = E1/kinEnergyProj;
dfsig *=dv/(-4.*std::log(y)/3.-4.*(1.-y)/3.+0.5*(1.-y*y));
dCrossEprod+=dfsig*cross1/dE;
}
*/
}
return dCrossEprod;
}
////////////////////////////////////////////////////////////////////////////////
//
G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond2(
const G4Material* aMaterial,
G4double kinEnergyProj, // kinetic energy of the primary particle before the interaction
G4double kinEnergyProd // kinetic energy of the secondary particle
)
{
G4double dCrossEprod=0.;
G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){
G4double dEdX1 = 0.0;
G4double dEdX2 = 0.0;
G4double E1=kinEnergyProd;
G4double E2=kinEnergyProd*1.001;
G4double dE=(E2-E1);
//G4double dum=0.;
dEdX1 = theDirectEMModel->ComputeDEDXPerVolume(aMaterial,G4Electron::Electron(),kinEnergyProj,E1);
dEdX2 = theDirectEMModel->ComputeDEDXPerVolume(aMaterial,G4Electron::Electron(),kinEnergyProj,E2);
dCrossEprod=(dEdX2-dEdX1)/dE/E1;
}
return dCrossEprod;
}
////////////////////////////////////////////////////////////////////////////////
//
G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond3(
const G4Material* aMaterial,
G4double kinEnergyProj, // kinetic energy of the primary particle before the interaction
G4double kinEnergyProd // kinetic energy of the secondary particle
)
{
return G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToSecond(aMaterial,
kinEnergyProj, // kinetic energy of the primary particle before the interaction
kinEnergyProd);
}
////////////////////////////////////////////////////////////////////////////////
//
G4double G4AdjointBremsstrahlungModel::SupressionFunction(const G4Material* material,
G4double kineticEnergy, G4double gammaEnergy)
{
// supression due to the LPM effect+polarisation of the medium/
// supression due to the polarisation alone
G4double totEnergy = kineticEnergy+electron_mass_c2 ;
G4double totEnergySquare = totEnergy*totEnergy ;
G4double LPMEnergy = LPMconstant*(material->GetRadlen()) ;
G4double gammaEnergySquare = gammaEnergy*gammaEnergy ;
G4double electronDensity = material->GetElectronDensity();
G4double sp = gammaEnergySquare/
(gammaEnergySquare+MigdalConstant*totEnergySquare*electronDensity);
G4double supr = 1.0;
if (theLPMflag) {
G4double s2lpm = LPMEnergy*gammaEnergy/totEnergySquare;
if (s2lpm < 1.) {
G4double LPMgEnergyLimit = totEnergySquare/LPMEnergy ;
G4double LPMgEnergyLimit2 = LPMgEnergyLimit*LPMgEnergyLimit;
G4double splim = LPMgEnergyLimit2/
(LPMgEnergyLimit2+MigdalConstant*totEnergySquare*electronDensity);
G4double w = 1.+1./splim ;
if ((1.-sp) < 1.e-6) w = s2lpm*(3.-sp);
else w = s2lpm*(1.+1./sp);
supr = (std::sqrt(w*w+4.*s2lpm)-w)/(std::sqrt(w*w+4.)-w) ;
supr /= sp;
}
}
return supr;
}
////////////////////////////////////////////////////////////////////////////////
@@ -491,59 +83,34 @@ G4double G4AdjointBremsstrahlungModel::SupressionFunction(const G4Material* mate
void G4AdjointBremsstrahlungModel::SampleSecondaries(const G4Track& aTrack,
G4bool IsScatProjToProjCase,
G4ParticleChange* fParticleChange)
{
{
if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange);
//G4cout<<"Adjoint Brem"<<std::endl;
const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
size_t ind=0;
if (UseMatrixPerElement ) { //Select Material
std::vector<double>* CS_Vs_Element = &CS_Vs_ElementForScatProjToProjCase;
if ( !IsScatProjToProjCase) CS_Vs_Element = &CS_Vs_ElementForProdToProjCase;
G4double rand_var= G4UniformRand();
G4double SumCS=0.;
for (size_t i=0;i<CS_Vs_Element->size();i++){
SumCS+=(*CS_Vs_Element)[i];
if (rand_var<=SumCS/lastCS){
ind=i;
break;
}
}
}
else {
ind = currentMaterialIndex;
}
const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
//Elastic inverse scattering modified compared to general G4VEmAdjointModel
//---------------------------
G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy();
//G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
return;
}
//Sample secondary energy
//-----------------------
G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(ind,
adjointPrimKinEnergy,
IsScatProjToProjCase);
G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy,
IsScatProjToProjCase);
//Weight correction
//-----------------------
CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(), adjointPrimKinEnergy,projectileKinEnergy);
CorrectPostStepWeight(fParticleChange,
aTrack.GetWeight(),
adjointPrimKinEnergy,
projectileKinEnergy,
IsScatProjToProjCase);
//Kinematic
//---------
G4double projectileM0 = electron_mass_c2;
G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
G4double projectileP = std::sqrt(projectileP2);
@@ -579,86 +146,249 @@ void G4AdjointBremsstrahlungModel::SampleSecondaries(const G4Track& aTrack,
if (!IsScatProjToProjCase && CorrectWeightMode){ //kill the primary and add a secondary
if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
fParticleChange->ProposeTrackStatus(fStopAndKill);
fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
//G4cout<<"projectileMomentum "<<projectileMomentum<<std::endl;
}
else {
fParticleChange->ProposeEnergy(projectileKinEnergy);
fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
//G4cout<<"projectileMomentum "<<projectileMomentum<<std::endl;
}
}
////////////////////////////////////////////////////////////////////////////////
//
void G4AdjointBremsstrahlungModel::DefineDirectBremModel(G4eBremsstrahlungModel* aModel)
{theDirectBremModel=aModel;
DefineDirectEMModel(aModel);
}
////////////////////////////////////////////////////////////////////////////////
//
void G4AdjointBremsstrahlungModel::InitialiseParameters()
{
static const G4double
ah10 = 4.67733E+00, ah11 =-6.19012E-01, ah12 = 2.02225E-02,
ah20 =-7.34101E+00, ah21 = 1.00462E+00, ah22 =-3.20985E-02,
ah30 = 2.93119E+00, ah31 =-4.03761E-01, ah32 = 1.25153E-02;
void G4AdjointBremsstrahlungModel::RapidSampleSecondaries(const G4Track& aTrack,
G4bool IsScatProjToProjCase,
G4ParticleChange* fParticleChange)
{
static const G4double
bh10 = 4.23071E+00, bh11 =-6.10995E-01, bh12 = 1.95531E-02,
bh20 =-7.12527E+00, bh21 = 9.69160E-01, bh22 =-2.74255E-02,
bh30 = 2.69925E+00, bh31 =-3.63283E-01, bh32 = 9.55316E-03;
/* static const G4double
al00 =-2.05398E+00, al01 = 2.38815E-02, al02 = 5.25483E-04,
al10 =-7.69748E-02, al11 =-6.91499E-02, al12 = 2.22453E-03,
al20 = 4.06463E-02, al21 =-1.01281E-02, al22 = 3.40919E-04;
static const G4double
bl00 = 1.04133E+00, bl01 =-9.43291E-03, bl02 =-4.54758E-04,
bl10 = 1.19253E-01, bl11 = 4.07467E-02, bl12 =-1.30718E-03,
bl20 =-1.59391E-02, bl21 = 7.27752E-03, bl22 =-1.94405E-04;*/
const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy();
if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
return;
}
const G4ElementTable* theElementTable = G4Element::GetElementTable();
FZ.clear();
ah1.clear();
ah2.clear();
ah3.clear();
bh1.clear();
bh2.clear();
bh3.clear();
al0.clear();
al1.clear();
al2.clear();
bl0.clear();
bl1.clear();
bl2.clear();
SigmaPerAtom.clear();
for (size_t j=0; j<theElementTable->size();j++){
G4double projectileKinEnergy =0.;
G4double gammaEnergy=0.;
G4double diffCSUsed=0.;
if (!IsScatProjToProjCase){
gammaEnergy=adjointPrimKinEnergy;
G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
G4double Emin= GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);;
if (Emin>=Emax) return;
projectileKinEnergy=Emin*std::pow(Emax/Emin,G4UniformRand());
diffCSUsed=lastCZ/projectileKinEnergy;
}
else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
G4double Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond);
if (Emin>=Emax) return;
G4double f1=(Emin-adjointPrimKinEnergy)/Emin;
G4double f2=(Emax-adjointPrimKinEnergy)/Emax/f1;
//G4cout<<"f1 and f2 "<<f1<<'\t'<<f2<<G4endl;
projectileKinEnergy=adjointPrimKinEnergy/(1.-f1*std::pow(f2,G4UniformRand()));
gammaEnergy=projectileKinEnergy-adjointPrimKinEnergy;
diffCSUsed=lastCZ*adjointPrimKinEnergy/projectileKinEnergy/gammaEnergy;
G4Element* anElement=(*theElementTable)[j];
G4double lnZ = 3.*(anElement->GetIonisation()->GetlogZ3());
FZ.push_back(lnZ* (4.- 0.55*lnZ));
G4double ZZ = anElement->GetIonisation()->GetZZ3();
ah1.push_back(ah10 + ZZ* (ah11 + ZZ* ah12));
ah2.push_back(ah20 + ZZ* (ah21 + ZZ* ah22));
ah3.push_back(ah30 + ZZ* (ah31 + ZZ* ah32));
}
//Weight correction
//-----------------------
//First w_corr is set to the ratio between adjoint total CS and fwd total CS
G4double w_corr=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
bh1.push_back(bh10 + ZZ* (bh11 + ZZ* bh12));
bh2.push_back(bh20 + ZZ* (bh21 + ZZ* bh22));
bh3.push_back(bh30 + ZZ* (bh31 + ZZ* bh32));
/*SigmaPerAtom.push_back(theDirectEMModel->ComputeCrossSectionPerAtom(
theDirectPrimaryPartDef,GetHighEnergyLimit()/2.,
anElement->GetZ(),1.,GetLowEnergyLimit(),1.e20));*/
//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
//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.
//Basically any other differential CS diffCS could be used here (example Penelope).
G4double diffCS = DiffCrossSectionPerVolumePrimToSecond(currentMaterial, projectileKinEnergy, gammaEnergy);
w_corr*=diffCS/diffCSUsed;
G4double new_weight = aTrack.GetWeight()*w_corr;
fParticleChange->SetParentWeightByProcess(false);
fParticleChange->SetSecondaryWeightByProcess(false);
fParticleChange->ProposeParentWeight(new_weight);
//Kinematic
//---------
G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
G4double projectileP = std::sqrt(projectileP2);
//Angle of the gamma direction with the projectile taken from G4eBremsstrahlungModel
//------------------------------------------------
G4double u;
const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
if (9./(9.+d) > G4UniformRand()) u = - std::log(G4UniformRand()*G4UniformRand())/a1;
else u = - std::log(G4UniformRand()*G4UniformRand())/a2;
G4double theta = u*electron_mass_c2/projectileTotalEnergy;
G4double sint = std::sin(theta);
G4double cost = std::cos(theta);
G4double phi = twopi * G4UniformRand() ;
G4ThreeVector projectileMomentum;
projectileMomentum=G4ThreeVector(std::cos(phi)*sint,std::sin(phi)*sint,cost)*projectileP; //gamma frame
if (IsScatProjToProjCase) {//the adjoint primary is the scattered e-
G4ThreeVector gammaMomentum = (projectileTotalEnergy-adjointPrimTotalEnergy)*G4ThreeVector(0.,0.,1.);
G4ThreeVector dirProd=projectileMomentum-gammaMomentum;
G4double cost1 = std::cos(dirProd.angle(projectileMomentum));
G4double sint1 = std::sqrt(1.-cost1*cost1);
projectileMomentum=G4ThreeVector(std::cos(phi)*sint1,std::sin(phi)*sint1,cost1)*projectileP;
}
projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection());
if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
fParticleChange->ProposeTrackStatus(fStopAndKill);
fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
}
else {
fParticleChange->ProposeEnergy(projectileKinEnergy);
fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
}
}
////////////////////////////////////////////////////////////////////////////////
//
G4AdjointBremsstrahlungModel::~G4AdjointBremsstrahlungModel()
{;}
////////////////////////////////////////////////////////////////////////////////
//
G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond(const G4Material* aMaterial,
G4double kinEnergyProj, // kinetic energy of the primary particle before the interaction
G4double kinEnergyProd // kinetic energy of the secondary particle
)
{/*if (UsePenelopeModel && !isPenelopeModelInitialised) {
theEmModelManagerForFwdModels->Initialise(G4Electron::Electron(),G4Gamma::Gamma(),1.,0);
isPenelopeModelInitialised =true;
}
*/
return DiffCrossSectionPerVolumePrimToSecondApproximated2(aMaterial,
kinEnergyProj,
kinEnergyProd);
/*return G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToSecond(aMaterial,
kinEnergyProj,
kinEnergyProd);*/
}
////////////////////////////////////////////////////////////////////////////////
//
G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecondApproximated1(
const G4Material* aMaterial,
G4double kinEnergyProj, // kinetic energy of the primary particle before the interaction
G4double kinEnergyProd // kinetic energy of the secondary particle
)
{
G4double dCrossEprod=0.;
G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
//In this approximation we consider that the secondary gammas are sampled with 1/Egamma energy distribution
//This is what is applied in the discrete standard model before the rejection test that make a cooerction
//The application of the same rejection function is not possble here.
//The differentiation of the CS over Ecut does not produce neither a good differential CS. That is due to the
// fact that in the discrete model the differential CS and the integrated CS are both fitted but separatly and
// therefore do not allow a correct numerical differentiation of the integrated CS to get the differential one.
// In the future we plan to use the brem secondary spectra from the G4Penelope implementation
if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){
G4double sigma=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,1.*keV);
dCrossEprod=sigma/kinEnergyProd/std::log(kinEnergyProj/keV);
}
return dCrossEprod;
}
////////////////////////////////////////////////////////////////////////////////
//
G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecondApproximated2(
const G4Material* material,
G4double kinEnergyProj, // kinetic energy of the primary particle before the interaction
G4double kinEnergyProd // kinetic energy of the secondary particle
)
{
//In this approximation we derive the direct cross section over Tcut=gamma energy, en after apply the Migdla correction factor
//used in the direct model
G4double dCrossEprod=0.;
const G4ElementVector* theElementVector = material->GetElementVector();
const double* theAtomNumDensityVector = material->GetAtomicNumDensityVector();
G4double dum=0.;
G4double E1=kinEnergyProd,E2=kinEnergyProd*1.001;
G4double dE=E2-E1;
for (size_t i=0; i<material->GetNumberOfElements(); i++) {
G4double C1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,(*theElementVector)[i]->GetZ(),dum ,E1);
G4double C2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,(*theElementVector)[i]->GetZ(),dum,E2);
dCrossEprod += theAtomNumDensityVector[i] * (C1-C2)/dE;
}
//Now the Migdal correction
G4double totalEnergy = kinEnergyProj+electron_mass_c2 ;
G4double kp2 = MigdalConstant*totalEnergy*totalEnergy
*(material->GetElectronDensity());
G4double MigdalFactor = 1./(1.+kp2/(kinEnergyProd*kinEnergyProd)); // its seems that the factor used in the CS compuation i the direct
//model is different than the one used in the secondary sampling by a
//factor (1.+kp2) To be checked!
dCrossEprod*=MigdalFactor;
return dCrossEprod;
}
////////////////////////////////////////////////////////////////////////////////
//
G4double G4AdjointBremsstrahlungModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
G4double primEnergy,
G4bool IsScatProjToProjCase)
{/* if (UsePenelopeModel && !isPenelopeModelInitialised) {
theEmModelManagerForFwdModels->Initialise(G4Electron::Electron(),G4Gamma::Gamma(),1.,0);
isPenelopeModelInitialised =true;
}
*/
if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase);
DefineCurrentMaterial(aCouple);
G4double Cross=0.;
lastCZ=theDirectEMModel->CrossSectionPerVolume(aCouple->GetMaterial(),theDirectPrimaryPartDef,100.*MeV,100.*MeV/std::exp(1.));//this give the constant above
if (!IsScatProjToProjCase ){
G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy);
G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy);
if (Emax_proj>Emin_proj && primEnergy > currentTcutForDirectSecond) Cross= lastCZ*std::log(Emax_proj/Emin_proj);
}
else {
G4double Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy);
G4double Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,currentTcutForDirectSecond);
if (Emax_proj>Emin_proj) Cross= lastCZ*std::log((Emax_proj-primEnergy)*Emin_proj/Emax_proj/(Emin_proj-primEnergy));
}
return Cross;
}
source/processes/electromagnetic/adjoint/src/G4AdjointCSManager.cc
+330 -209
View File
@@ -23,6 +23,9 @@
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
// $Id: G4AdjointCSManager.cc,v 1.5 2009/11/20 10:31:20 ldesorgh Exp $
// GEANT4 tag $Name: geant4-09-03 $
//
#include "G4AdjointCSManager.hh"
#include "G4AdjointCSMatrix.hh"
#include "G4AdjointInterpolator.hh"
@@ -39,10 +42,14 @@
#include "G4PhysicsTableHelper.hh"
#include "G4Electron.hh"
#include "G4Gamma.hh"
#include "G4Proton.hh"
#include "G4AdjointElectron.hh"
#include "G4AdjointGamma.hh"
#include "G4AdjointProton.hh"
#include "G4ProductionCutsTable.hh"
#include "G4ProductionCutsTable.hh"
#include <fstream>
#include <iomanip>
G4AdjointCSManager* G4AdjointCSManager::theInstance = 0;
@@ -64,18 +71,32 @@ G4AdjointCSManager::G4AdjointCSManager()
theTotalAdjointSigmaTableVector.clear();
listOfForwardEmProcess.clear();
listOfForwardEnergyLossProcess.clear();
theListOfAdjointParticlesInAction.clear();
theListOfAdjointParticlesInAction.clear();
EminForFwdSigmaTables.clear();
EminForAdjSigmaTables.clear();
EkinofFwdSigmaMax.clear();
EkinofAdjSigmaMax.clear();
Tmin=0.1*keV;
Tmax=100.*TeV;
nbins=240;
nbins=360; //probably this should be decrease, that was choosen to avoid error in the CS value closed to CS jump.(For example at Tcut)
RegisterAdjointParticle(G4AdjointElectron::AdjointElectron());
RegisterAdjointParticle(G4AdjointGamma::AdjointGamma());
RegisterAdjointParticle(G4AdjointProton::AdjointProton());
verbose = 1;
lastPartDefForCS =0;
LastEkinForCS =0;
LastCSCorrectionFactor =1.;
consider_continuous_weight_correction =true;
consider_poststep_weight_correction =false;
forward_CS_mode = true;
currentParticleDef = 0;
theAdjIon = 0;
theFwdIon = 0;
}
///////////////////////////////////////////////////////
@@ -95,7 +116,7 @@ void G4AdjointCSManager::RegisterEmProcess(G4VEmProcess* aProcess, G4ParticleDef
G4ParticleDefinition* anAdjPartDef = GetAdjointParticleEquivalent(aFwdPartDef);
if (anAdjPartDef && aProcess){
RegisterAdjointParticle(anAdjPartDef);
int index=-1;
G4int index=-1;
for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
if (anAdjPartDef->GetParticleName() == theListOfAdjointParticlesInAction[i]->GetParticleName()) index=i;
@@ -110,7 +131,7 @@ void G4AdjointCSManager::RegisterEnergyLossProcess(G4VEnergyLossProcess* aProces
G4ParticleDefinition* anAdjPartDef = GetAdjointParticleEquivalent(aFwdPartDef);
if (anAdjPartDef && aProcess){
RegisterAdjointParticle(anAdjPartDef);
int index=-1;
G4int index=-1;
for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
if (anAdjPartDef->GetParticleName() == theListOfAdjointParticlesInAction[i]->GetParticleName()) index=i;
}
@@ -120,7 +141,7 @@ void G4AdjointCSManager::RegisterEnergyLossProcess(G4VEnergyLossProcess* aProces
///////////////////////////////////////////////////////
//
void G4AdjointCSManager::RegisterAdjointParticle(G4ParticleDefinition* aPartDef)
{ int index=-1;
{ G4int index=-1;
for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
if (aPartDef->GetParticleName() == theListOfAdjointParticlesInAction[i]->GetParticleName()) index=i;
}
@@ -131,6 +152,11 @@ void G4AdjointCSManager::RegisterAdjointParticle(G4ParticleDefinition* aPartDef)
theTotalAdjointSigmaTableVector.push_back(new G4PhysicsTable);
listOfForwardEmProcess.push_back(new std::vector<G4VEmProcess*>());
theListOfAdjointParticlesInAction.push_back(aPartDef);
EminForFwdSigmaTables.push_back(std::vector<G4double> ());
EminForAdjSigmaTables.push_back(std::vector<G4double> ());
EkinofFwdSigmaMax.push_back(std::vector<G4double> ());
EkinofAdjSigmaMax.push_back(std::vector<G4double> ());
}
}
///////////////////////////////////////////////////////
@@ -161,9 +187,11 @@ void G4AdjointCSManager::BuildCrossSectionMatrices()
theAdjointCSMatricesForProdToProj.clear();
const G4ElementTable* theElementTable = G4Element::GetElementTable();
const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
G4cout<<"========== Computation of cross section matrices for adjoint models =========="<<G4endl;
for (size_t i=0; i<listOfAdjointEMModel.size();i++){
G4VEmAdjointModel* aModel =listOfAdjointEMModel[i];
G4cout<<"Build adjoint cross section matrices for "<<aModel->GetName()<<std::endl;
G4cout<<"Build adjoint cross section matrices for "<<aModel->GetName()<<G4endl;
if (aModel->GetUseMatrix()){
std::vector<G4AdjointCSMatrix*>* aListOfMat1 = new std::vector<G4AdjointCSMatrix*>();
std::vector<G4AdjointCSMatrix*>* aListOfMat2 = new std::vector<G4AdjointCSMatrix*>();
@@ -172,17 +200,17 @@ void G4AdjointCSManager::BuildCrossSectionMatrices()
if (aModel->GetUseMatrixPerElement()){
if (aModel->GetUseOnlyOneMatrixForAllElements()){
std::vector<G4AdjointCSMatrix*>
two_matrices=BuildCrossSectionsMatricesForAGivenModelAndElement(aModel,1, 1, 10);
two_matrices=BuildCrossSectionsMatricesForAGivenModelAndElement(aModel,1, 1, 80);
aListOfMat1->push_back(two_matrices[0]);
aListOfMat2->push_back(two_matrices[1]);
}
else {
for (size_t j=0; j<theElementTable->size();j++){
G4Element* anElement=(*theElementTable)[j];
G4int Z = G4int(anElement->GetZ());
G4int A = G4int(anElement->GetA());
G4int Z = int(anElement->GetZ());
G4int A = int(anElement->GetA());
std::vector<G4AdjointCSMatrix*>
two_matrices=BuildCrossSectionsMatricesForAGivenModelAndElement(aModel,Z, A, 10);
two_matrices=BuildCrossSectionsMatricesForAGivenModelAndElement(aModel,Z, A, 40);
aListOfMat1->push_back(two_matrices[0]);
aListOfMat2->push_back(two_matrices[1]);
}
@@ -192,7 +220,7 @@ void G4AdjointCSManager::BuildCrossSectionMatrices()
for (size_t j=0; j<theMaterialTable->size();j++){
G4Material* aMaterial=(*theMaterialTable)[j];
std::vector<G4AdjointCSMatrix*>
two_matrices=BuildCrossSectionsMatricesForAGivenModelAndMaterial(aModel,aMaterial, 10);
two_matrices=BuildCrossSectionsMatricesForAGivenModelAndMaterial(aModel,aMaterial, 40);
aListOfMat1->push_back(two_matrices[0]);
aListOfMat2->push_back(two_matrices[1]);
}
@@ -202,14 +230,19 @@ void G4AdjointCSManager::BuildCrossSectionMatrices()
theAdjointCSMatricesForScatProjToProj.push_back(*aListOfMat2);
aModel->SetCSMatrices(aListOfMat1, aListOfMat2);
}
else { std::vector<G4AdjointCSMatrix*> two_empty_matrices;
else { G4cout<<"The model "<<aModel->GetName()<<" does not use cross section matrices"<<G4endl;
std::vector<G4AdjointCSMatrix*> two_empty_matrices;
theAdjointCSMatricesForProdToProj.push_back(two_empty_matrices);
theAdjointCSMatricesForScatProjToProj.push_back(two_empty_matrices);
}
}
G4cout<<"All adjoint cross section matrices are built "<<std::endl;
G4cout<<" All adjoint cross section matrices are computed!"<<G4endl;
G4cout<<"======================================================================"<<G4endl;
CrossSectionMatrixesAreBuilt = true;
}
@@ -220,36 +253,106 @@ void G4AdjointCSManager::BuildTotalSigmaTables()
const G4ProductionCutsTable* theCoupleTable= G4ProductionCutsTable::GetProductionCutsTable();
for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
G4ParticleDefinition* thePartDef = theListOfAdjointParticlesInAction[i];
DefineCurrentParticle(thePartDef);
theTotalForwardSigmaTableVector[i]->clearAndDestroy();
theTotalAdjointSigmaTableVector[i]->clearAndDestroy();
EminForFwdSigmaTables[i].clear();
EminForAdjSigmaTables[i].clear();
EkinofFwdSigmaMax[i].clear();
EkinofAdjSigmaMax[i].clear();
//G4cout<<thePartDef->GetParticleName();
for (size_t j=0;j<theCoupleTable->GetTableSize();j++){
const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
/*
G4String file_name1=couple->GetMaterial()->GetName()+"_"+thePartDef->GetParticleName()+"_adj_totCS.txt";
G4String file_name2=couple->GetMaterial()->GetName()+"_"+thePartDef->GetParticleName()+"_fwd_totCS.txt";
std::fstream FileOutputAdjCS(file_name1, std::ios::out);
std::fstream FileOutputFwdCS(file_name2, std::ios::out);
FileOutputAdjCS<<std::setiosflags(std::ios::scientific);
FileOutputAdjCS<<std::setprecision(6);
FileOutputFwdCS<<std::setiosflags(std::ios::scientific);
FileOutputFwdCS<<std::setprecision(6);
*/
//make first the total fwd CS table for FwdProcess
G4PhysicsVector* aVector = new G4PhysicsLogVector(Tmin, Tmax, nbins);
G4bool Emin_found=false;
size_t ind=0;
G4double sigma_max =0.;
G4double e_sigma_max =0.;
for(size_t l=0; l<aVector->GetVectorLength(); l++) {
G4double totCS=0;
G4double totCS=0.;
G4double e=aVector->GetLowEdgeEnergy(l);
for (size_t k=0; k<listOfForwardEmProcess[i]->size(); k++){
totCS+=(*listOfForwardEmProcess[i])[k]->GetLambda(e, couple);
}
for (size_t k=0; k<listOfForwardEnergyLossProcess[i]->size(); k++){
totCS+=(*listOfForwardEnergyLossProcess[i])[k]->GetLambda(e, couple);
if (thePartDef == theAdjIon) { // e is considered already as the scaled energy
size_t mat_index = couple->GetIndex();
G4VEmModel* currentModel = (*listOfForwardEnergyLossProcess[i])[k]->SelectModelForMaterial(e,mat_index);
G4double chargeSqRatio = currentModel->GetChargeSquareRatio(theFwdIon,couple->GetMaterial(),e/massRatio);
(*listOfForwardEnergyLossProcess[i])[k]->SetDynamicMassCharge(massRatio,chargeSqRatio);
}
G4double e1=e/massRatio;
totCS+=(*listOfForwardEnergyLossProcess[i])[k]->GetLambda(e1, couple);
}
//G4cout<<totCS<<std::endl;
aVector->PutValue(l,totCS);
if (totCS>sigma_max){
sigma_max=totCS;
e_sigma_max = e;
}
//FileOutputFwdCS<<e<<'\t'<<totCS<<G4endl;
if (totCS>0 && !Emin_found) {
EminForFwdSigmaTables[i].push_back(e);
Emin_found=true;
}
}
//FileOutputFwdCS.close();
EkinofFwdSigmaMax[i].push_back(e_sigma_max);
if(!Emin_found) EminForFwdSigmaTables[i].push_back(Tmax);
theTotalForwardSigmaTableVector[i]->push_back(aVector);
Emin_found=false;
sigma_max=0;
e_sigma_max =0.;
ind=0;
G4PhysicsVector* aVector1 = new G4PhysicsLogVector(Tmin, Tmax, nbins);
for(size_t l=0; l<aVector->GetVectorLength(); l++) {
G4double e=aVector->GetLowEdgeEnergy(l);
G4double totCS =ComputeTotalAdjointCS(couple,thePartDef,e);
//G4cout<<totCS<<std::endl;
G4double totCS =ComputeTotalAdjointCS(couple,thePartDef,e*0.9999999/massRatio); //massRatio needed for ions
aVector1->PutValue(l,totCS);
if (totCS>sigma_max){
sigma_max=totCS;
e_sigma_max = e;
}
//FileOutputAdjCS<<e<<'\t'<<totCS<<G4endl;
if (totCS>0 && !Emin_found) {
EminForAdjSigmaTables[i].push_back(e);
Emin_found=true;
}
}
}
//FileOutputAdjCS.close();
EkinofAdjSigmaMax[i].push_back(e_sigma_max);
if(!Emin_found) EminForAdjSigmaTables[i].push_back(Tmax);
theTotalAdjointSigmaTableVector[i]->push_back(aVector1);
}
@@ -261,14 +364,9 @@ void G4AdjointCSManager::BuildTotalSigmaTables()
G4double G4AdjointCSManager::GetTotalAdjointCS(G4ParticleDefinition* aPartDef, G4double Ekin,
const G4MaterialCutsCouple* aCouple)
{ DefineCurrentMaterial(aCouple);
int index=-1;
for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
if (aPartDef == theListOfAdjointParticlesInAction[i]) index=i;
}
if (index == -1) return 0.;
DefineCurrentParticle(aPartDef);
G4bool b;
return (((*theTotalAdjointSigmaTableVector[index])[currentMatIndex])->GetValue(Ekin, b));
return (((*theTotalAdjointSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(Ekin*massRatio, b));
@@ -278,70 +376,142 @@ G4double G4AdjointCSManager::GetTotalAdjointCS(G4ParticleDefinition* aPartDef, G
G4double G4AdjointCSManager::GetTotalForwardCS(G4ParticleDefinition* aPartDef, G4double Ekin,
const G4MaterialCutsCouple* aCouple)
{ DefineCurrentMaterial(aCouple);
int index=-1;
for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
if (aPartDef == theListOfAdjointParticlesInAction[i]) index=i;
}
if (index == -1) return 0.;
DefineCurrentParticle(aPartDef);
G4bool b;
return (((*theTotalForwardSigmaTableVector[index])[currentMatIndex])->GetValue(Ekin, b));
return (((*theTotalForwardSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(Ekin*massRatio, b));
}
///////////////////////////////////////////////////////
//
void G4AdjointCSManager::GetEminForTotalCS(G4ParticleDefinition* aPartDef,
const G4MaterialCutsCouple* aCouple, G4double& emin_adj, G4double& emin_fwd)
{ DefineCurrentMaterial(aCouple);
DefineCurrentParticle(aPartDef);
emin_adj = EminForAdjSigmaTables[currentParticleIndex][currentMatIndex]/massRatio;
emin_fwd = EminForFwdSigmaTables[currentParticleIndex][currentMatIndex]/massRatio;
}
///////////////////////////////////////////////////////
//
void G4AdjointCSManager::GetMaxFwdTotalCS(G4ParticleDefinition* aPartDef,
const G4MaterialCutsCouple* aCouple, G4double& e_sigma_max, G4double& sigma_max)
{ DefineCurrentMaterial(aCouple);
DefineCurrentParticle(aPartDef);
e_sigma_max = EkinofFwdSigmaMax[currentParticleIndex][currentMatIndex];
G4bool b;
sigma_max =((*theTotalForwardSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(e_sigma_max, b);
e_sigma_max/=massRatio;
}
///////////////////////////////////////////////////////
//
void G4AdjointCSManager::GetMaxAdjTotalCS(G4ParticleDefinition* aPartDef,
const G4MaterialCutsCouple* aCouple, G4double& e_sigma_max, G4double& sigma_max)
{ DefineCurrentMaterial(aCouple);
DefineCurrentParticle(aPartDef);
e_sigma_max = EkinofAdjSigmaMax[currentParticleIndex][currentMatIndex];
G4bool b;
sigma_max =((*theTotalAdjointSigmaTableVector[currentParticleIndex])[currentMatIndex])->GetValue(e_sigma_max, b);
e_sigma_max/=massRatio;
}
///////////////////////////////////////////////////////
//
G4double G4AdjointCSManager::GetCrossSectionCorrection(G4ParticleDefinition* aPartDef,G4double PreStepEkin,const G4MaterialCutsCouple* aCouple, G4bool& fwd_is_used,
G4double& fwd_TotCS)
{ G4double corr_fac = 1.;
if (forward_CS_mode) {
fwd_TotCS=PrefwdCS;
if (LastEkinForCS != PreStepEkin || aPartDef != lastPartDefForCS || aCouple!=currentCouple) {
DefineCurrentMaterial(aCouple);
PreadjCS = GetTotalAdjointCS(aPartDef, PreStepEkin,aCouple);
PrefwdCS = GetTotalForwardCS(aPartDef, PreStepEkin,aCouple);
LastEkinForCS = PreStepEkin;
lastPartDefForCS = aPartDef;
if (PrefwdCS >0. && PreadjCS >0.) {
forward_CS_is_used = true;
LastCSCorrectionFactor = PrefwdCS/PreadjCS;
}
else {
forward_CS_is_used = false;
LastCSCorrectionFactor = 1.;
}
}
corr_fac =LastCSCorrectionFactor;
}
else {
forward_CS_is_used = false;
LastCSCorrectionFactor = 1.;
}
fwd_TotCS=PrefwdCS;
fwd_is_used = forward_CS_is_used;
return corr_fac;
}
///////////////////////////////////////////////////////
//
G4double G4AdjointCSManager::GetContinuousWeightCorrection(G4ParticleDefinition* aPartDef, G4double PreStepEkin,G4double AfterStepEkin,
const G4MaterialCutsCouple* aCouple, G4double step_length)
{ //G4double fwdCS = GetTotalForwardCS(aPartDef, AfterStepEkin,aCouple);
G4double corr_fac = 1.;
if (consider_continuous_weight_correction) {
G4double adjCS = GetTotalAdjointCS(aPartDef, PreStepEkin,aCouple);
G4double PrefwdCS;
PrefwdCS = GetTotalForwardCS(aPartDef, PreStepEkin,aCouple);
G4double fwdCS = GetTotalForwardCS(aPartDef, (AfterStepEkin+PreStepEkin)/2.,aCouple);
G4cout<<adjCS<<'\t'<<fwdCS<<std::endl;
//if (aPartDef ==G4AdjointGamma::AdjointGamma()) G4cout<<adjCS<<'\t'<<fwdCS<<std::endl;
/*if (adjCS >0 ) corr_fac = std::exp((PrefwdCS-fwdCS)*step_length);
else corr_fac = std::exp(-fwdCS*step_length);*/
corr_fac *=std::exp((adjCS-fwdCS)*step_length);
corr_fac=std::max(corr_fac,1.e-6);
corr_fac *=PreStepEkin/AfterStepEkin;
{ G4double corr_fac = 1.;
//return corr_fac;
//G4double after_adjCS = GetTotalAdjointCS(aPartDef, AfterStepEkin,aCouple);
G4double after_fwdCS = GetTotalForwardCS(aPartDef, AfterStepEkin,aCouple);
G4double pre_adjCS = GetTotalAdjointCS(aPartDef, PreStepEkin,aCouple);
if (!forward_CS_is_used || pre_adjCS ==0. || after_fwdCS==0.) {
forward_CS_is_used=false;
G4double pre_fwdCS = GetTotalForwardCS(aPartDef, PreStepEkin,aCouple);
corr_fac *=std::exp((pre_adjCS-pre_fwdCS)*step_length);
LastCSCorrectionFactor = 1.;
}
G4cout<<"Cont "<<corr_fac<<std::endl;
G4cout<<"Ekin0 "<<PreStepEkin<<std::endl;
G4cout<<"Ekin1 "<<AfterStepEkin<<std::endl;
G4cout<<"step_length "<<step_length<<std::endl;
else {
LastCSCorrectionFactor = after_fwdCS/pre_adjCS;
}
return corr_fac;
}
///////////////////////////////////////////////////////
//
G4double G4AdjointCSManager::GetPostStepWeightCorrection(G4ParticleDefinition* , G4ParticleDefinition* ,
G4double ,G4double ,
const G4MaterialCutsCouple* )
{ G4double corr_fac = 1.;
if (consider_poststep_weight_correction) {
/*G4double fwdCS = GetTotalForwardCS(aSecondPartDef, EkinPrim,aCouple);
G4double adjCS = GetTotalAdjointCS(aPrimPartDef, EkinPrim,aCouple);*/
//G4double fwd1CS = GetTotalForwardCS(aPrimPartDef, EkinPrim,aCouple);
//if (adjCS>0 && fwd1CS>0) adjCS = fwd1CS;
//corr_fac =fwdCS*EkinSecond/adjCS/EkinPrim;
//corr_fac = adjCS/fwdCS;
}
return corr_fac;
G4double G4AdjointCSManager::GetPostStepWeightCorrection( )
{//return 1.;
return 1./LastCSCorrectionFactor;
}
///////////////////////////////////////////////////////
//
double G4AdjointCSManager::ComputeAdjointCS(G4Material* aMaterial,
G4double G4AdjointCSManager::ComputeAdjointCS(G4Material* aMaterial,
G4VEmAdjointModel* aModel,
G4double PrimEnergy,
G4double Tcut,
G4bool IsScatProjToProjCase,
std::vector<double>& CS_Vs_Element)
std::vector<G4double>& CS_Vs_Element)
{
G4double EminSec=0;
G4double EmaxSec=0;
if (IsScatProjToProjCase){
EminSec= aModel->GetSecondAdjEnergyMinForScatProjToProjCase(PrimEnergy,Tcut);
EmaxSec= aModel->GetSecondAdjEnergyMaxForScatProjToProjCase(PrimEnergy);
}
else if (PrimEnergy > Tcut || !aModel->GetApplyCutInRange()) {
EminSec= aModel->GetSecondAdjEnergyMinForProdToProjCase(PrimEnergy);
EmaxSec= aModel->GetSecondAdjEnergyMaxForProdToProjCase(PrimEnergy);
}
if (EminSec >= EmaxSec) return 0.;
G4bool need_to_compute=false;
if ( aMaterial!= lastMaterial || PrimEnergy != lastPrimaryEnergy || Tcut != lastTcut){
lastMaterial =aMaterial;
@@ -380,7 +550,7 @@ double G4AdjointCSManager::ComputeAdjointCS(G4Material* aMaterial,
listOfIsScatProjToProjCase.push_back(IsScatProjToProjCase);
CS_Vs_Element.clear();
if (!aModel->GetUseMatrix()){
return aModel->AdjointCrossSection(currentCouple,PrimEnergy,IsScatProjToProjCase);
CS_Vs_Element.push_back(aModel->AdjointCrossSection(currentCouple,PrimEnergy,IsScatProjToProjCase));
}
@@ -396,18 +566,9 @@ double G4AdjointCSManager::ComputeAdjointCS(G4Material* aMaterial,
if (PrimEnergy > Tlow)
CS = ComputeAdjointCS(PrimEnergy,theCSMatrix,Tlow);
G4double factor=0.;
for (size_t i=0;i<n_el;i++){
size_t ind_el = aMaterial->GetElement(i)->GetIndex();
for (size_t i=0;i<n_el;i++){ //this could be computed only once
//size_t ind_el = aMaterial->GetElement(i)->GetIndex();
factor+=aMaterial->GetElement(i)->GetZ()*aMaterial->GetVecNbOfAtomsPerVolume()[i];
G4AdjointCSMatrix* theCSMatrix;
if (IsScatProjToProjCase){
theCSMatrix=theAdjointCSMatricesForScatProjToProj[ind_model][ind_el];
}
else theCSMatrix=theAdjointCSMatricesForProdToProj[ind_model][ind_el];
//G4double CS =0.;
//G4cout<<CS<<std::endl;
}
CS *=factor;
CS_Vs_Element.push_back(CS);
@@ -416,7 +577,7 @@ double G4AdjointCSManager::ComputeAdjointCS(G4Material* aMaterial,
else {
for (size_t i=0;i<n_el;i++){
size_t ind_el = aMaterial->GetElement(i)->GetIndex();
//G4cout<<aMaterial->GetName()<<std::endl;
//G4cout<<aMaterial->GetName()<<G4endl;
G4AdjointCSMatrix* theCSMatrix;
if (IsScatProjToProjCase){
theCSMatrix=theAdjointCSMatricesForScatProjToProj[ind_model][ind_el];
@@ -425,7 +586,7 @@ double G4AdjointCSManager::ComputeAdjointCS(G4Material* aMaterial,
G4double CS =0.;
if (PrimEnergy > Tlow)
CS = ComputeAdjointCS(PrimEnergy,theCSMatrix,Tlow);
//G4cout<<CS<<std::endl;
//G4cout<<CS<<G4endl;
CS_Vs_Element.push_back(CS*(aMaterial->GetVecNbOfAtomsPerVolume()[i]));
}
}
@@ -452,18 +613,10 @@ double G4AdjointCSManager::ComputeAdjointCS(G4Material* aMaterial,
G4double CS=0;
for (size_t i=0;i<CS_Vs_Element.size();i++){
CS+=CS_Vs_Element[i];
}
return CS;
CS+=CS_Vs_Element[i]; //We could put the progressive sum of the CS instead of the CS of an element itself
}
return CS;
}
///////////////////////////////////////////////////////
//
@@ -472,7 +625,7 @@ G4Element* G4AdjointCSManager::SampleElementFromCSMatrices(G4Material* aMaterial
G4double PrimEnergy,
G4double Tcut,
G4bool IsScatProjToProjCase)
{ std::vector<double> CS_Vs_Element;
{ std::vector<G4double> CS_Vs_Element;
G4double CS = ComputeAdjointCS(aMaterial,aModel,PrimEnergy,Tcut,IsScatProjToProjCase,CS_Vs_Element);
G4double rand_var= G4UniformRand();
G4double SumCS=0.;
@@ -497,59 +650,40 @@ G4double G4AdjointCSManager::ComputeTotalAdjointCS(const G4MaterialCutsCouple* a
G4double Ekin)
{
G4double TotalCS=0.;
// G4ParticleDefinition* theDirPartDef = GetForwardParticleEquivalent(aPartDef);
DefineCurrentMaterial(aCouple);
/* size_t idx=-1;
if (theDirPartDef->GetParticleName() == "gamma") idx = 0;
else if (theDirPartDef->GetParticleName() == "e-") idx = 1;
else if (theDirPartDef->GetParticleName() == "e+") idx = 2;
//THe tCut computation is wrong this should be on Tcut per model the secondary determioming the Tcut
const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
//G4cout<<aVec<<std::endl;
G4double Tcut =(*aVec)[aCouple->GetIndex()];*/
//G4cout<<"Tcut "<<Tcut<<std::endl;
//G4cout<<(*aVec)[0]<<std::endl;
// G4double Tcut =converters[idx]->Convert(Rcut,aCouple->GetMaterial());
std::vector<double> CS_Vs_Element;
std::vector<G4double> CS_Vs_Element;
for (size_t i=0; i<listOfAdjointEMModel.size();i++){
/*G4ParticleDefinition* theDirSecondPartDef =
GetForwardParticleEquivalent(listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectSecondaryParticleDefinition());
*/
G4double Tlow=0;
if (!listOfAdjointEMModel[i]->GetApplyCutInRange()) Tlow =listOfAdjointEMModel[i]->GetLowEnergyLimit();
else {
G4ParticleDefinition* theDirSecondPartDef =
GetForwardParticleEquivalent(listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectSecondaryParticleDefinition());
G4int idx=-1;
size_t idx=56;
if (theDirSecondPartDef->GetParticleName() == "gamma") idx = 0;
else if (theDirSecondPartDef->GetParticleName() == "e-") idx = 1;
else if (theDirSecondPartDef->GetParticleName() == "e+") idx = 2;
const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
Tlow =(*aVec)[aCouple->GetIndex()];
if (idx <56) {
const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
Tlow =(*aVec)[aCouple->GetIndex()];
}
}
if ( Ekin<=listOfAdjointEMModel[i]->GetHighEnergyLimit() && Ekin>=listOfAdjointEMModel[i]->GetLowEnergyLimit()){
if (aPartDef == listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectPrimaryParticleDefinition()){
//G4cout<<"Yes1 before "<<std::endl;
TotalCS += ComputeAdjointCS(currentMaterial,
listOfAdjointEMModel[i],
Ekin, Tlow,true,CS_Vs_Element);
//G4cout<<"Yes1 "<<Ekin<<'\t'<<TotalCS<<std::endl;
Ekin, Tlow,true,CS_Vs_Element);
}
if (aPartDef == listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectSecondaryParticleDefinition()){
TotalCS += ComputeAdjointCS(currentMaterial,
listOfAdjointEMModel[i],
Ekin, Tlow,false, CS_Vs_Element);
//G4cout<<"Yes2 "<<TotalCS<<std::endl;
}
}
@@ -562,7 +696,7 @@ G4double G4AdjointCSManager::ComputeTotalAdjointCS(const G4MaterialCutsCouple* a
//
std::vector<G4AdjointCSMatrix*>
G4AdjointCSManager::BuildCrossSectionsMatricesForAGivenModelAndElement(G4VEmAdjointModel* aModel,G4int Z,G4int A,
int nbin_pro_decade)
G4int nbin_pro_decade)
{
G4AdjointCSMatrix* theCSMatForProdToProjBackwardScattering = new G4AdjointCSMatrix(false);
G4AdjointCSMatrix* theCSMatForScatProjToProjBackwardScattering = new G4AdjointCSMatrix(true);
@@ -576,36 +710,25 @@ G4AdjointCSManager::BuildCrossSectionsMatricesForAGivenModelAndElement(G4VEmAdjo
if (aModel->GetSecondPartOfSameType() )EkinMaxForProd =EkinMaxForProd/2.;
//Product to projectile backward scattering
//-----------------------------------------
G4double fE=std::pow(10.,1./nbin_pro_decade);
G4double E2=std::pow(10.,G4double( G4int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
G4double E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
G4double E1=EkinMin;
while (E1 <EkinMaxForProd){
E1=std::max(EkinMin,E2);
E1=std::min(EkinMaxForProd,E1);
std::vector< std::vector< G4double >* > aMat= aModel->ComputeAdjointCrossSectionVectorPerAtomForSecond(E1,Z,A,nbin_pro_decade);
std::vector< std::vector< double>* > aMat= aModel->ComputeAdjointCrossSectionVectorPerAtomForSecond(E1,Z,A,nbin_pro_decade);
if (aMat.size()>=2) {
std::vector< G4double >* log_ESecVec=aMat[0];
std::vector< G4double >* log_CSVec=aMat[1];
std::vector< double>* log_ESecVec=aMat[0];
std::vector< double>* log_CSVec=aMat[1];
G4double log_adjointCS=log_CSVec->back();
//normalise CSVec such that it becomes a probability vector
/*for (size_t j=0;j<log_CSVec->size();j++) (*log_CSVec)[j]=(*log_CSVec)[j]-log_adjointCS;
(*log_CSVec)[0]=-90.;*/
for (size_t j=0;j<log_CSVec->size();j++) {
//G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<std::endl;
if (j==0) (*log_CSVec)[j] = 0.;
else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS));
//G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<std::endl;
for (size_t j=0;j<log_CSVec->size();j++) {
if (j==0) (*log_CSVec)[j] = 0.;
else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS) +1e-50);
}
(*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-1.;
(*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
theCSMatForProdToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
}
E1=E2;
@@ -615,24 +738,22 @@ G4AdjointCSManager::BuildCrossSectionsMatricesForAGivenModelAndElement(G4VEmAdjo
//Scattered projectile to projectile backward scattering
//-----------------------------------------
E2=std::pow(10.,G4double( G4int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
E1=EkinMin;
while (E1 <EkinMaxForScat){
E1=std::max(EkinMin,E2);
E1=std::min(EkinMaxForScat,E1);
std::vector< std::vector< G4double >* > aMat= aModel->ComputeAdjointCrossSectionVectorPerAtomForScatProj(E1,Z,A,nbin_pro_decade);
std::vector< std::vector< double>* > aMat= aModel->ComputeAdjointCrossSectionVectorPerAtomForScatProj(E1,Z,A,nbin_pro_decade);
if (aMat.size()>=2) {
std::vector< G4double >* log_ESecVec=aMat[0];
std::vector< G4double >* log_CSVec=aMat[1];
std::vector< double>* log_ESecVec=aMat[0];
std::vector< double>* log_CSVec=aMat[1];
G4double log_adjointCS=log_CSVec->back();
//normalise CSVec such that it becomes a probability vector
for (size_t j=0;j<log_CSVec->size();j++) {
//G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<std::endl;
if (j==0) (*log_CSVec)[j] = 0.;
else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS));
//G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<std::endl;
if (j==0) (*log_CSVec)[j] = 0.;
else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS)+1e-50);
}
(*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-1.;
(*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
theCSMatForScatProjToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
}
E1=E2;
@@ -640,19 +761,13 @@ G4AdjointCSManager::BuildCrossSectionsMatricesForAGivenModelAndElement(G4VEmAdjo
}
std::vector<G4AdjointCSMatrix*> res;
res.clear();
res.push_back(theCSMatForProdToProjBackwardScattering);
res.push_back(theCSMatForScatProjToProjBackwardScattering);
#ifdef TEST_MODE
/*
G4String file_name;
std::stringstream astream;
G4String str_Z;
@@ -661,14 +776,8 @@ G4AdjointCSManager::BuildCrossSectionsMatricesForAGivenModelAndElement(G4VEmAdjo
theCSMatForProdToProjBackwardScattering->Write(aModel->GetName()+G4String("_CSMat_Z")+str_Z+"_ProdToProj.txt");
theCSMatForScatProjToProjBackwardScattering->Write(aModel->GetName()+G4String("_CSMat_Z")+str_Z+"_ScatProjToProj.txt");
/*G4AdjointCSMatrix* aMat1 = new G4AdjointCSMatrix(false);
G4AdjointCSMatrix* aMat2 = new G4AdjointCSMatrix(true);
aMat1->Read(G4String("test_Z")+str_Z+"_1.txt");
aMat2->Read(G4String("test_Z")+str_Z+"_2.txt");
aMat1->Write(G4String("test_Z")+str_Z+"_11.txt");
aMat2->Write(G4String("test_Z")+str_Z+"_22.txt"); */
#endif
*/
return res;
@@ -701,25 +810,25 @@ G4AdjointCSManager::BuildCrossSectionsMatricesForAGivenModelAndMaterial(G4VEmAdj
//Product to projectile backward scattering
//-----------------------------------------
G4double fE=std::pow(10.,1./nbin_pro_decade);
G4double E2=std::pow(10.,G4double( G4int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
G4double E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
G4double E1=EkinMin;
while (E1 <EkinMaxForProd){
E1=std::max(EkinMin,E2);
E1=std::min(EkinMaxForProd,E1);
std::vector< std::vector< G4double >* > aMat= aModel->ComputeAdjointCrossSectionVectorPerVolumeForSecond(aMaterial,E1,nbin_pro_decade);
std::vector< std::vector< double>* > aMat= aModel->ComputeAdjointCrossSectionVectorPerVolumeForSecond(aMaterial,E1,nbin_pro_decade);
if (aMat.size()>=2) {
std::vector< G4double >* log_ESecVec=aMat[0];
std::vector< G4double >* log_CSVec=aMat[1];
std::vector< double>* log_ESecVec=aMat[0];
std::vector< double>* log_CSVec=aMat[1];
G4double log_adjointCS=log_CSVec->back();
//normalise CSVec such that it becomes a probability vector
for (size_t j=0;j<log_CSVec->size();j++) {
//G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<std::endl;
//G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<G4endl;
if (j==0) (*log_CSVec)[j] = 0.;
else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS));
//G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<std::endl;
//G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<G4endl;
}
(*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-1.;
(*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
theCSMatForProdToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
}
@@ -732,24 +841,25 @@ G4AdjointCSManager::BuildCrossSectionsMatricesForAGivenModelAndMaterial(G4VEmAdj
//Scattered projectile to projectile backward scattering
//-----------------------------------------
E2=std::pow(10.,G4double( G4int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
E2=std::pow(10.,double( int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
E1=EkinMin;
while (E1 <EkinMaxForScat){
E1=std::max(EkinMin,E2);
E1=std::min(EkinMaxForScat,E1);
std::vector< std::vector< G4double >* > aMat= aModel->ComputeAdjointCrossSectionVectorPerVolumeForScatProj(aMaterial,E1,nbin_pro_decade);
std::vector< std::vector< double>* > aMat= aModel->ComputeAdjointCrossSectionVectorPerVolumeForScatProj(aMaterial,E1,nbin_pro_decade);
if (aMat.size()>=2) {
std::vector< G4double >* log_ESecVec=aMat[0];
std::vector< G4double >* log_CSVec=aMat[1];
std::vector< double>* log_ESecVec=aMat[0];
std::vector< double>* log_CSVec=aMat[1];
G4double log_adjointCS=log_CSVec->back();
for (size_t j=0;j<log_CSVec->size();j++) {
//G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<std::endl;
//G4cout<<"CSMan1 "<<(*log_CSVec)[j]<<G4endl;
if (j==0) (*log_CSVec)[j] = 0.;
else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS));
//G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<std::endl;
//G4cout<<"CSMan2 "<<(*log_CSVec)[j]<<G4endl;if (theAdjPartDef->GetParticleName() == "adj_gamma") return G4Gamma::Gamma();
}
(*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-1.;
(*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
theCSMatForScatProjToProjBackwardScattering->AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
}
@@ -769,10 +879,10 @@ G4AdjointCSManager::BuildCrossSectionsMatricesForAGivenModelAndMaterial(G4VEmAdj
res.push_back(theCSMatForProdToProjBackwardScattering);
res.push_back(theCSMatForScatProjToProjBackwardScattering);
#ifdef TEST_MODE
/*
theCSMatForProdToProjBackwardScattering->Write(aModel->GetName()+"_CSMat_"+aMaterial->GetName()+"_ProdToProj.txt");
theCSMatForScatProjToProjBackwardScattering->Write(aModel->GetName()+"_CSMat_"+aMaterial->GetName()+"_ScatProjToProj.txt");
#endif
*/
return res;
@@ -785,7 +895,10 @@ G4AdjointCSManager::BuildCrossSectionsMatricesForAGivenModelAndMaterial(G4VEmAdj
G4ParticleDefinition* G4AdjointCSManager::GetAdjointParticleEquivalent(G4ParticleDefinition* theFwdPartDef)
{
if (theFwdPartDef->GetParticleName() == "e-") return G4AdjointElectron::AdjointElectron();
if (theFwdPartDef->GetParticleName() == "gamma") return G4AdjointGamma::AdjointGamma();
else if (theFwdPartDef->GetParticleName() == "gamma") return G4AdjointGamma::AdjointGamma();
else if (theFwdPartDef->GetParticleName() == "proton") return G4AdjointProton::AdjointProton();
else if (theFwdPartDef ==theFwdIon) return theAdjIon;
return 0;
}
///////////////////////////////////////////////////////
@@ -793,7 +906,9 @@ G4ParticleDefinition* G4AdjointCSManager::GetAdjointParticleEquivalent(G4Particl
G4ParticleDefinition* G4AdjointCSManager::GetForwardParticleEquivalent(G4ParticleDefinition* theAdjPartDef)
{
if (theAdjPartDef->GetParticleName() == "adj_e-") return G4Electron::Electron();
if (theAdjPartDef->GetParticleName() == "adj_gamma") return G4Gamma::Gamma();
else if (theAdjPartDef->GetParticleName() == "adj_gamma") return G4Gamma::Gamma();
else if (theAdjPartDef->GetParticleName() == "adj_proton") return G4Proton::Proton();
else if (theAdjPartDef == theAdjIon) return theFwdIon;
return 0;
}
///////////////////////////////////////////////////////
@@ -804,25 +919,43 @@ void G4AdjointCSManager::DefineCurrentMaterial(const G4MaterialCutsCouple* coupl
currentCouple = const_cast<G4MaterialCutsCouple*> (couple);
currentMaterial = const_cast<G4Material*> (couple->GetMaterial());
currentMatIndex = couple->GetIndex();
//G4cout<<"Index material "<<currentMatIndex<<std::endl;
lastPartDefForCS =0;
LastEkinForCS =0;
LastCSCorrectionFactor =1.;
}
}
///////////////////////////////////////////////////////
//
void G4AdjointCSManager::DefineCurrentParticle(const G4ParticleDefinition* aPartDef)
{
if(aPartDef != currentParticleDef) {
currentParticleDef= const_cast< G4ParticleDefinition* > (aPartDef);
massRatio=1;
if (aPartDef == theAdjIon) massRatio = proton_mass_c2/aPartDef->GetPDGMass();
currentParticleIndex=1000000;
for (size_t i=0;i<theListOfAdjointParticlesInAction.size();i++){
if (aPartDef == theListOfAdjointParticlesInAction[i]) currentParticleIndex=i;
}
}
}
///////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////
//
double G4AdjointCSManager::ComputeAdjointCS(G4double aPrimEnergy,G4AdjointCSMatrix*
G4double G4AdjointCSManager::ComputeAdjointCS(G4double aPrimEnergy,G4AdjointCSMatrix*
anAdjointCSMatrix,G4double Tcut)
{
std::vector< G4double > *theLogPrimEnergyVector = anAdjointCSMatrix->GetLogPrimEnergyVector();
std::vector< double> *theLogPrimEnergyVector = anAdjointCSMatrix->GetLogPrimEnergyVector();
if (theLogPrimEnergyVector->size() ==0){
G4cout<<"No data are contained in the given AdjointCSMatrix!"<<std::endl;
G4cout<<"The sampling procedure will be stopped."<<std::endl;
G4cout<<"No data are contained in the given AdjointCSMatrix!"<<G4endl;
G4cout<<"The s"<<G4endl;
return 0.;
}
//G4cout<<"A prim/Tcut "<<aPrimEnergy<<'\t'<<Tcut<<std::endl;
G4double log_Tcut = std::log(Tcut);
G4double log_E =std::log(aPrimEnergy);
@@ -833,35 +966,23 @@ double G4AdjointCSManager::ComputeAdjointCS(G4double aPrimEnergy,G4AdjointCSMatr
G4AdjointInterpolator* theInterpolator=G4AdjointInterpolator::GetInstance();
size_t ind =theInterpolator->FindPositionForLogVector(log_E,*theLogPrimEnergyVector);
//G4cout<<"Prim energy "<<(*thePrimEnergyVector)[0]<<std::endl;
//G4cout<<"Prim energy[ind]"<<(*thePrimEnergyVector)[ind]<<std::endl;
//G4cout<<"Prim energy ind"<<ind<<std::endl;
G4double aLogPrimEnergy1,aLogPrimEnergy2;
G4double aLogCS1,aLogCS2;
G4double log01,log02;
std::vector< G4double>* aLogSecondEnergyVector1 =0;
std::vector< G4double>* aLogSecondEnergyVector2 =0;
std::vector< G4double>* aLogProbVector1=0;
std::vector< G4double>* aLogProbVector2=0;
std::vector< double>* aLogSecondEnergyVector1 =0;
std::vector< double>* aLogSecondEnergyVector2 =0;
std::vector< double>* aLogProbVector1=0;
std::vector< double>* aLogProbVector2=0;
std::vector< size_t>* aLogProbVectorIndex1=0;
std::vector< size_t>* aLogProbVectorIndex2=0;
anAdjointCSMatrix->GetData(ind, aLogPrimEnergy1,aLogCS1,log01, aLogSecondEnergyVector1,aLogProbVector1,aLogProbVectorIndex1);
anAdjointCSMatrix->GetData(ind+1, aLogPrimEnergy2,aLogCS2,log02, aLogSecondEnergyVector2,aLogProbVector2,aLogProbVectorIndex2);
//G4cout<<"aSecondEnergyVector1.size() "<<aSecondEnergyVector1->size()<<std::endl;
//G4cout<<aSecondEnergyVector1<<std::endl;
//G4cout<<"aSecondEnergyVector2.size() "<<aSecondEnergyVector2->size()<<std::endl;
if (anAdjointCSMatrix->IsScatProjToProjCase()){ //case where the Tcut plays a role
G4double log_minimum_prob1, log_minimum_prob2;
//G4cout<<aSecondEnergyVector1->size()<<std::endl;
log_minimum_prob1=theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector1,*aLogProbVector1);
log_minimum_prob2=theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector2,*aLogProbVector2);
//G4cout<<"minimum_prob1 "<< std::exp(log_minimum_prob1)<<std::endl;
//G4cout<<"minimum_prob2 "<< std::exp(log_minimum_prob2)<<std::endl;
//G4cout<<"Tcut "<<std::endl;
aLogCS1+= log_minimum_prob1;
aLogCS2+= log_minimum_prob2;
}
source/processes/electromagnetic/adjoint/src/G4AdjointCSMatrix.cc
+23 -28
View File
@@ -23,12 +23,14 @@
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
// $Id: G4AdjointCSMatrix.cc,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
// GEANT4 tag $Name: geant4-09-03 $
//
#include "G4AdjointCSMatrix.hh"
#include <iomanip>
#include <fstream>
#include "G4AdjointInterpolator.hh"
///////////////////////////////////////////////////////
//
G4AdjointCSMatrix::G4AdjointCSMatrix(G4bool aBool){
@@ -63,27 +65,26 @@ void G4AdjointCSMatrix::Clear()
}
///////////////////////////////////////////////////////
//
void G4AdjointCSMatrix::AddData(G4double aLogPrimEnergy,G4double aLogCS, std::vector< G4double>* aLogSecondEnergyVector,
std::vector< G4double>* aLogProbVector,size_t n_pro_decade){
void G4AdjointCSMatrix::AddData(G4double aLogPrimEnergy,G4double aLogCS, std::vector< double>* aLogSecondEnergyVector,
std::vector< double>* aLogProbVector,size_t n_pro_decade){
G4AdjointInterpolator* theInterpolator=G4AdjointInterpolator::GetInstance();
//Add this time we consider that the energy are given monotically
//At this time we consider that the energy is increasing monotically
theLogPrimEnergyVector.push_back(aLogPrimEnergy);
theLogCrossSectionVector.push_back(aLogCS);
theLogSecondEnergyMatrix.push_back(aLogSecondEnergyVector);
//G4cout<<"Test Add Data "<<this<<'\t'<<aSecondEnergyVector->size()<<std::endl;
//G4cout<<theSecondEnergyMatrix.size()<<std::endl;
theLogProbMatrix.push_back(aLogProbVector);
//G4cout<<"Test Add Data 1 "<<this<<'\t'<<aSecondEnergyVector->size()<<std::endl;
//G4cout<<theSecondEnergyMatrix.size()<<std::endl;
std::vector< size_t>* aLogProbVectorIndex = 0;
dlog =0;
if (n_pro_decade > 0 && aLogProbVector->size()>0) {
aLogProbVectorIndex = new std::vector< size_t>();
dlog=std::log(10.)/n_pro_decade;
G4double log_val = int(std::min((*aLogProbVector)[0],aLogProbVector->back())/dlog)*dlog;
log0Vector.push_back(log_val);
while(log_val<0.) {
aLogProbVectorIndex->push_back(theInterpolator->FindPosition(log_val,(*aLogProbVector)));
log_val+=dlog;
@@ -101,22 +102,16 @@ void G4AdjointCSMatrix::Clear()
}
///////////////////////////////////////////////////////
//
bool G4AdjointCSMatrix::GetData(unsigned int i, G4double& aLogPrimEnergy,G4double& aLogCS,G4double& log0, std::vector< G4double>*& aLogSecondEnergyVector,
std::vector< G4double>*& aLogProbVector, std::vector< size_t>*& aLogProbVectorIndex)
G4bool G4AdjointCSMatrix::GetData(unsigned int i, G4double& aLogPrimEnergy,G4double& aLogCS,G4double& log0, std::vector< double>*& aLogSecondEnergyVector,
std::vector< double>*& aLogProbVector, std::vector< size_t>*& aLogProbVectorIndex)
{ if (i>= nb_of_PrimEnergy) return false;
//G4cout<<"Test Get Data "<<std::endl;
//G4cout<<"Test Get Data "<<G4endl;
aLogPrimEnergy = theLogPrimEnergyVector[i];
aLogCS = theLogCrossSectionVector[i];
aLogSecondEnergyVector = theLogSecondEnergyMatrix[i];
//G4cout<<"Test Get Data "<<this<<'\t'<<theSecondEnergyMatrix[i]->size()<<std::endl;
//G4cout<<"Test Get Data "<<this<<'\t'<<aSecondEnergyVector->size()<<std::endl;
//G4cout<<"Test Get Data "<<this<<'\t'<<aSecondEnergyVector<<std::endl;
aLogProbVector = theLogProbMatrix[i];
aLogProbVectorIndex = theLogProbMatrixIndex[i];
log0=log0Vector[i];
//G4cout<<"Test Get Data 1 "<<this<<'\t'<<theProbMatrix[i]->size()<<std::endl;
//G4cout<<"Test Get Data 1 "<<this<<'\t'<<aProbVector->size()<<std::endl;
//G4cout<<"Test Get Data 1 "<<this<<'\t'<<aLogProbVectorIndex<<std::endl;
return true;
}
@@ -126,33 +121,33 @@ void G4AdjointCSMatrix::Write(G4String file_name)
{ std::fstream FileOutput(file_name, std::ios::out);
FileOutput<<std::setiosflags(std::ios::scientific);
FileOutput<<std::setprecision(6);
FileOutput<<theLogPrimEnergyVector.size()<<std::endl;
FileOutput<<theLogPrimEnergyVector.size()<<G4endl;
for (size_t i=0;i<theLogPrimEnergyVector.size();i++){
FileOutput<<std::exp(theLogPrimEnergyVector[i])/MeV<<'\t'<<std::exp(theLogCrossSectionVector[i])<<std::endl;
FileOutput<<std::exp(theLogPrimEnergyVector[i])/MeV<<'\t'<<std::exp(theLogCrossSectionVector[i])<<G4endl;
size_t j1=0;
FileOutput<<theLogSecondEnergyMatrix[i]->size()<<std::endl;
FileOutput<<theLogSecondEnergyMatrix[i]->size()<<G4endl;
for (size_t j=0;j<theLogSecondEnergyMatrix[i]->size();j++){
FileOutput<<std::exp((*theLogSecondEnergyMatrix[i])[j]);
j1++;
if (j1<10) FileOutput<<'\t';
else {
FileOutput<<std::endl;
FileOutput<<G4endl;
j1=0;
}
}
if (j1>0) FileOutput<<std::endl;
if (j1>0) FileOutput<<G4endl;
j1=0;
FileOutput<<theLogProbMatrix[i]->size()<<std::endl;
FileOutput<<theLogProbMatrix[i]->size()<<G4endl;
for (size_t j=0;j<theLogProbMatrix[i]->size();j++){
FileOutput<<std::exp((*theLogProbMatrix[i])[j]);
j1++;
if (j1<10) FileOutput<<'\t';
else {
FileOutput<<std::endl;
FileOutput<<G4endl;
j1=0;
}
}
if (j1>0) FileOutput<<std::endl;
if (j1>0) FileOutput<<G4endl;
}
@@ -176,8 +171,8 @@ void G4AdjointCSMatrix::Read(G4String file_name)
theLogPrimEnergyVector.push_back(E);
theLogCrossSectionVector.push_back(CS);
FileOutput>>n2;
theLogSecondEnergyMatrix.push_back(new std::vector<double>());
theLogProbMatrix.push_back(new std::vector<double>());
theLogSecondEnergyMatrix.push_back(new std::vector<G4double>());
theLogProbMatrix.push_back(new std::vector<G4double>());
for (size_t j=0; j<n2;j++){
G4double E1;
source/processes/electromagnetic/adjoint/src/G4AdjointComptonModel.cc
+191 -35
View File
@@ -23,6 +23,9 @@
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
// $Id: G4AdjointComptonModel.cc,v 1.6 2009/12/16 17:50:03 gunter Exp $
// GEANT4 tag $Name: geant4-09-03 $
//
#include "G4AdjointComptonModel.hh"
#include "G4AdjointCSManager.hh"
@@ -33,6 +36,7 @@
#include "G4AdjointElectron.hh"
#include "G4AdjointGamma.hh"
#include "G4Gamma.hh"
#include "G4KleinNishinaCompton.hh"
////////////////////////////////////////////////////////////////////////////////
@@ -43,12 +47,12 @@ G4AdjointComptonModel::G4AdjointComptonModel():
{ SetApplyCutInRange(false);
SetUseMatrix(true);
SetUseMatrixPerElement(true);
SetIsIonisation(false);
SetUseOnlyOneMatrixForAllElements(true);
theAdjEquivOfDirectPrimPartDef =G4AdjointGamma::AdjointGamma();
theAdjEquivOfDirectSecondPartDef=G4AdjointElectron::AdjointElectron();
theDirectPrimaryPartDef=G4Gamma::Gamma();
second_part_of_same_type=false;
theDirectEMModel=new G4KleinNishinaCompton(G4Gamma::Gamma(),"ComptonDirectModel");
}
////////////////////////////////////////////////////////////////////////////////
@@ -61,7 +65,7 @@ void G4AdjointComptonModel::SampleSecondaries(const G4Track& aTrack,
G4bool IsScatProjToProjCase,
G4ParticleChange* fParticleChange)
{
if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange);
//A recall of the compton scattering law is
//Egamma2=Egamma1/(1+(Egamma1/E0_electron)(1.-cos_th))
@@ -70,27 +74,16 @@ void G4AdjointComptonModel::SampleSecondaries(const G4Track& aTrack,
const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
//DefineCurrentMaterial(aTrack->GetMaterialCutsCouple());
size_t ind= 0;
//Elastic inverse scattering //not correct in all the cases
//---------------------------------------------------------
G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
//G4cout<<adjointPrimKinEnergy<<std::endl;
if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
return;
}
}
//Sample secondary energy
//-----------------------
G4double gammaE1;
gammaE1 = SampleAdjSecEnergyFromCSMatrix(ind,
adjointPrimKinEnergy,
gammaE1 = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy,
IsScatProjToProjCase);
@@ -107,7 +100,7 @@ void G4AdjointComptonModel::SampleSecondaries(const G4Track& aTrack,
//Cos th
//-------
// G4cout<<"Compton scattering "<<gammaE1<<'\t'<<gammaE2<<std::endl;
// G4cout<<"Compton scattering "<<gammaE1<<'\t'<<gammaE2<<G4endl;
G4double cos_th = 1.+ electron_mass_c2*(1./gammaE1 -1./gammaE2);
if (!IsScatProjToProjCase) {
G4double p_elec=theAdjointPrimary->GetTotalMomentum();
@@ -115,7 +108,7 @@ void G4AdjointComptonModel::SampleSecondaries(const G4Track& aTrack,
}
G4double sin_th = 0.;
if (std::abs(cos_th)>1){
//G4cout<<"Problem in compton scattering with cos_th "<<cos_th<<std::endl;
//G4cout<<"Problem in compton scattering with cos_th "<<cos_th<<G4endl;
if (cos_th>0) {
cos_th=1.;
}
@@ -135,17 +128,21 @@ void G4AdjointComptonModel::SampleSecondaries(const G4Track& aTrack,
G4double phi =G4UniformRand()*2.*3.1415926;
G4ThreeVector gammaMomentum1 = gammaE1*G4ThreeVector(std::cos(phi)*sin_th,std::sin(phi)*sin_th,cos_th);
gammaMomentum1.rotateUz(dir_parallel);
// G4cout<<gamma0Energy<<'\t'<<gamma0Momentum<<std::endl;
// G4cout<<gamma0Energy<<'\t'<<gamma0Momentum<<G4endl;
//It is important to correct the weight of particles before adding the secondary
//------------------------------------------------------------------------------
CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(), adjointPrimKinEnergy,gammaE1);
CorrectPostStepWeight(fParticleChange,
aTrack.GetWeight(),
adjointPrimKinEnergy,
gammaE1,
IsScatProjToProjCase);
if (!IsScatProjToProjCase && CorrectWeightMode){ //kill the primary and add a secondary
if (!IsScatProjToProjCase){ //kill the primary and add a secondary
fParticleChange->ProposeTrackStatus(fStopAndKill);
fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,gammaMomentum1));
//G4cout<<"gamma0Momentum "<<gamma0Momentum<<std::endl;
//G4cout<<"gamma0Momentum "<<gamma0Momentum<<G4endl;
}
else {
fParticleChange->ProposeEnergy(gammaE1);
@@ -153,7 +150,127 @@ void G4AdjointComptonModel::SampleSecondaries(const G4Track& aTrack,
}
}
}
////////////////////////////////////////////////////////////////////////////////
//
void G4AdjointComptonModel::RapidSampleSecondaries(const G4Track& aTrack,
G4bool IsScatProjToProjCase,
G4ParticleChange* fParticleChange)
{
const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
return;
}
G4double diffCSUsed=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2;
G4double gammaE1=0.;
G4double gammaE2=0.;
if (!IsScatProjToProjCase){
G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
G4double Emin= GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);;
if (Emin>=Emax) return;
G4double f1=(Emin-adjointPrimKinEnergy)/Emin;
G4double f2=(Emax-adjointPrimKinEnergy)/Emax/f1;
gammaE1=adjointPrimKinEnergy/(1.-f1*std::pow(f2,G4UniformRand()));;
gammaE2=gammaE1-adjointPrimKinEnergy;
diffCSUsed= diffCSUsed*(1.+2.*std::log(1.+electron_mass_c2/adjointPrimKinEnergy))*adjointPrimKinEnergy/gammaE1/gammaE2;
}
else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
G4double Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond);
if (Emin>=Emax) return;
gammaE2 =adjointPrimKinEnergy;
gammaE1=Emin*std::pow(Emax/Emin,G4UniformRand());
diffCSUsed= diffCSUsed/gammaE1;
}
//Weight correction
//-----------------------
//First w_corr is set to the ratio between adjoint total CS and fwd total CS
G4double w_corr=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
//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
G4double diffCS = DiffCrossSectionPerAtomPrimToScatPrim(gammaE1, gammaE2,1,0.);
if (diffCS >0) diffCS /=G4direct_CS; // here we have the normalised diffCS
diffCS*=theDirectEMProcess->GetLambda(gammaE1,currentCouple);
//diffCS*=theDirectEMModel->CrossSectionPerVolume(currentMaterial,G4Gamma::Gamma(),gammaE1,0.,2.*gammaE1);
//G4cout<<"diffCS/diffCSUsed "<<diffCS/diffCSUsed<<'\t'<<gammaE1<<'\t'<<gammaE2<<G4endl;
w_corr*=diffCS/diffCSUsed;
G4double new_weight = aTrack.GetWeight()*w_corr;
fParticleChange->SetParentWeightByProcess(false);
fParticleChange->SetSecondaryWeightByProcess(false);
fParticleChange->ProposeParentWeight(new_weight);
//Cos th
//-------
G4double cos_th = 1.+ electron_mass_c2*(1./gammaE1 -1./gammaE2);
if (!IsScatProjToProjCase) {
G4double p_elec=theAdjointPrimary->GetTotalMomentum();
cos_th = (gammaE1 - gammaE2*cos_th)/p_elec;
}
G4double sin_th = 0.;
if (std::abs(cos_th)>1){
//G4cout<<"Problem in compton scattering with cos_th "<<cos_th<<G4endl;
if (cos_th>0) {
cos_th=1.;
}
else cos_th=-1.;
sin_th=0.;
}
else sin_th = std::sqrt(1.-cos_th*cos_th);
//gamma0 momentum
//--------------------
G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
G4double phi =G4UniformRand()*2.*3.1415926;
G4ThreeVector gammaMomentum1 = gammaE1*G4ThreeVector(std::cos(phi)*sin_th,std::sin(phi)*sin_th,cos_th);
gammaMomentum1.rotateUz(dir_parallel);
if (!IsScatProjToProjCase){ //kill the primary and add a secondary
fParticleChange->ProposeTrackStatus(fStopAndKill);
fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,gammaMomentum1));
//G4cout<<"gamma0Momentum "<<gamma0Momentum<<G4endl;
}
else {
fParticleChange->ProposeEnergy(gammaE1);
fParticleChange->ProposeMomentumDirection(gammaMomentum1.unit());
}
}
////////////////////////////////////////////////////////////////////////////////
//
//The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
@@ -178,9 +295,11 @@ G4double G4AdjointComptonModel::DiffCrossSectionPerAtomPrimToScatPrim(
G4double Z,
G4double )
{ //Based on Klein Nishina formula
//* In the forward case (see G4KleinNishinaModel) the cross section is parametrised while the secondaries are sampled from the
// In the forward case (see G4KleinNishinaModel) the cross section is parametrised while
// the secondaries are sampled from the
// Klein Nishida differential cross section
// The used diffrential cross section here is therefore the cross section multiplied by the normalidsed differential Klein Nishida cross section
// The used diffrential cross section here is therefore the cross section multiplied by the normalised
//differential Klein Nishida cross section
//Klein Nishida Cross Section
@@ -190,10 +309,10 @@ G4double G4AdjointComptonModel::DiffCrossSectionPerAtomPrimToScatPrim(
G4double gamEnergy1_max = gamEnergy0;
G4double gamEnergy1_min = gamEnergy0/one_plus_two_epsi;
if (gamEnergy1 >gamEnergy1_max || gamEnergy1<gamEnergy1_min) {
/*G4cout<<"the differential CS is null"<<std::endl;
G4cout<<gamEnergy0<<std::endl;
G4cout<<gamEnergy1<<std::endl;
G4cout<<gamEnergy1_min<<std::endl;*/
/*G4cout<<"the differential CS is null"<<G4endl;
G4cout<<gamEnergy0<<G4endl;
G4cout<<gamEnergy1<<G4endl;
G4cout<<gamEnergy1_min<<G4endl;*/
return 0.;
}
@@ -221,19 +340,20 @@ G4double G4AdjointComptonModel::DiffCrossSectionPerAtomPrimToScatPrim(
//Normalised to the CS used in G4
//-------------------------------
G4double G4direct_CS = theDirectEMModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(),
G4direct_CS = theDirectEMModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(),
gamEnergy0,
Z, 0., 0.,0.);
dCS_dE1 *= G4direct_CS/CS;
/* G4cout<<"the differential CS is not null"<<std::endl;
G4cout<<gamEnergy0<<std::endl;
G4cout<<gamEnergy1<<std::endl;*/
/* G4cout<<"the differential CS is not null"<<G4endl;
G4cout<<gamEnergy0<<G4endl;
G4cout<<gamEnergy1<<G4endl;*/
return dCS_dE1;
}
////////////////////////////////////////////////////////////////////////////////
//
G4double G4AdjointComptonModel::GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
@@ -250,3 +370,39 @@ G4double G4AdjointComptonModel::GetSecondAdjEnergyMinForProdToProjCase(G4double
G4double emin=half_e+term;
return emin;
}
////////////////////////////////////////////////////////////////////////////////
//
G4double G4AdjointComptonModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
G4double primEnergy,
G4bool IsScatProjToProjCase)
{
if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase);
DefineCurrentMaterial(aCouple);
G4double Cross=0.;
G4double Emax_proj =0.;
G4double Emin_proj =0.;
if (!IsScatProjToProjCase ){
Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy);
Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy);
if (Emax_proj>Emin_proj ){
Cross= std::log((Emax_proj-primEnergy)*Emin_proj/Emax_proj/(Emin_proj-primEnergy))
*(1.+2.*std::log(1.+electron_mass_c2/primEnergy));
}
}
else {
Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy);
Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,0.);
if (Emax_proj>Emin_proj) {
Cross = std::log(Emax_proj/Emin_proj);
//+0.5*primEnergy*primEnergy(1./(Emin_proj*Emin_proj) - 1./(Emax_proj*Emax_proj)); neglected at the moment
}
}
Cross*=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2;
lastCS=Cross;
return Cross;
}
source/processes/electromagnetic/adjoint/src/G4AdjointInterpolator.cc
+25 -22
View File
@@ -23,6 +23,9 @@
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
// $Id: G4AdjointInterpolator.cc,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
// GEANT4 tag $Name: geant4-09-03 $
//
#include "G4AdjointCSMatrix.hh"
#include "G4AdjointInterpolator.hh"
@@ -60,7 +63,7 @@ G4AdjointInterpolator::~G4AdjointInterpolator()
//
G4double G4AdjointInterpolator::LinearInterpolation(G4double& x,G4double& x1,G4double& x2,G4double& y1,G4double& y2)
{ G4double res = y1+ (x-x1)*(y2-y1)/(x2-x1);
//G4cout<<"Linear "<<res<<std::endl;
//G4cout<<"Linear "<<res<<G4endl;
return res;
}
///////////////////////////////////////////////////////
@@ -68,10 +71,10 @@ G4double G4AdjointInterpolator::LinearInterpolation(G4double& x,G4double& x1,G4d
G4double G4AdjointInterpolator::LogarithmicInterpolation(G4double& x,G4double& x1,G4double& x2,G4double& y1,G4double& y2)
{ if (y1<=0 || y2<=0 || x1<=0) return LinearInterpolation(x,x1,x2,y1,y2);
G4double B=std::log(y2/y1)/std::log(x2/x1);
//G4cout<<"x1,x2,y1,y2 "<<x1<<'\t'<<x2<<'\t'<<y1<<'\t'<<y2<<'\t'<<std::endl;
//G4cout<<"x1,x2,y1,y2 "<<x1<<'\t'<<x2<<'\t'<<y1<<'\t'<<y2<<'\t'<<G4endl;
G4double A=y1/std::pow(x1,B);
G4double res=A*std::pow(x,B);
// G4cout<<"Log "<<res<<std::endl;
// G4cout<<"Log "<<res<<G4endl;
return res;
}
///////////////////////////////////////////////////////
@@ -97,15 +100,15 @@ G4double G4AdjointInterpolator::Interpolation(G4double& x,G4double& x1,G4double&
return ExponentialInterpolation(x,x1,x2,y1,y2);
}
else {
//G4cout<<"The interpolation method that you invoked does not exist!"<<std::endl;
//G4cout<<"The interpolation method that you invoked does not exist!"<<G4endl;
return -1111111111.;
}
}
///////////////////////////////////////////////////////
//
size_t G4AdjointInterpolator::FindPosition(G4double& x,std::vector<double>& x_vec,size_t , size_t ) //only valid if x_vec is monotically increasing
size_t G4AdjointInterpolator::FindPosition(G4double& x,std::vector<G4double>& x_vec,size_t , size_t ) //only valid if x_vec is monotically increasing
{ //most rapid nethod could be used probably
//It is important to put std::vector<double>& such that the vector itself is used and not a copy
//It is important to put std::vector<G4double>& such that the vector itself is used and not a copy
size_t ndim = x_vec.size();
@@ -148,10 +151,10 @@ size_t G4AdjointInterpolator::FindPosition(G4double& x,std::vector<double>& x_v
///////////////////////////////////////////////////////
//
size_t G4AdjointInterpolator::FindPositionForLogVector(G4double& log_x,std::vector<double>& log_x_vec) //only valid if x_vec is monotically increasing
size_t G4AdjointInterpolator::FindPositionForLogVector(G4double& log_x,std::vector<G4double>& log_x_vec) //only valid if x_vec is monotically increasing
{ //most rapid nethod could be used probably
//It is important to put std::vector<double>& such that the vector itself is used and not a copy
//It is important to put std::vector<G4double>& such that the vector itself is used and not a copy
return FindPosition(log_x, log_x_vec);
if (log_x_vec.size()>3){
size_t ind=0;
G4double log_x1=log_x_vec[1];
@@ -169,17 +172,17 @@ size_t G4AdjointInterpolator::FindPositionForLogVector(G4double& log_x,std::vec
}
///////////////////////////////////////////////////////
//
G4double G4AdjointInterpolator::Interpolate(G4double& x,std::vector<double>& x_vec,std::vector<double>& y_vec,G4String InterPolMethod)
G4double G4AdjointInterpolator::Interpolate(G4double& x,std::vector<G4double>& x_vec,std::vector<G4double>& y_vec,G4String InterPolMethod)
{ size_t i=FindPosition(x,x_vec);
//G4cout<<i<<std::endl;
//G4cout<<x<<std::endl;
//G4cout<<x_vec[i]<<std::endl;
//G4cout<<i<<G4endl;
//G4cout<<x<<G4endl;
//G4cout<<x_vec[i]<<G4endl;
return Interpolation( x,x_vec[i],x_vec[i+1],y_vec[i],y_vec[i+1],InterPolMethod);
}
///////////////////////////////////////////////////////
//
G4double G4AdjointInterpolator::InterpolateWithIndexVector(G4double& x,std::vector<double>& x_vec,std::vector<double>& y_vec,
G4double G4AdjointInterpolator::InterpolateWithIndexVector(G4double& x,std::vector<G4double>& x_vec,std::vector<G4double>& y_vec,
std::vector<size_t>& index_vec,G4double x0, G4double dx) //only linear interpolation possible
{ size_t ind=0;
if (x>x0) ind=int((x-x0)/dx);
@@ -201,16 +204,16 @@ G4double G4AdjointInterpolator::InterpolateWithIndexVector(G4double& x,std::vect
///////////////////////////////////////////////////////
//
G4double G4AdjointInterpolator::InterpolateForLogVector(G4double& log_x,std::vector<double>& log_x_vec,std::vector<double>& log_y_vec)
G4double G4AdjointInterpolator::InterpolateForLogVector(G4double& log_x,std::vector<G4double>& log_x_vec,std::vector<G4double>& log_y_vec)
{ //size_t i=0;
size_t i=FindPositionForLogVector(log_x,log_x_vec);
/*G4cout<<"In interpolate "<<std::endl;
G4cout<<i<<std::endl;
G4cout<<log_x<<std::endl;
G4cout<<log_x_vec[i]<<std::endl;
G4cout<<log_x_vec[i+1]<<std::endl;
G4cout<<log_y_vec[i]<<std::endl;
G4cout<<log_y_vec[i+1]<<std::endl;*/
/*G4cout<<"In interpolate "<<G4endl;
G4cout<<i<<G4endl;
G4cout<<log_x<<G4endl;
G4cout<<log_x_vec[i]<<G4endl;
G4cout<<log_x_vec[i+1]<<G4endl;
G4cout<<log_y_vec[i]<<G4endl;
G4cout<<log_y_vec[i+1]<<G4endl;*/
G4double log_y=LinearInterpolation(log_x,log_x_vec[i],log_x_vec[i+1],log_y_vec[i],log_y_vec[i+1]);
return log_y;
source/processes/electromagnetic/adjoint/src/G4AdjointIonIonisationModel.cc
+363
View File
@@ -0,0 +1,363 @@
//
// ********************************************************************
// * License and Disclaimer *
// * *
// * The Geant4 software is copyright of the Copyright Holders of *
// * the Geant4 Collaboration. It is provided under the terms and *
// * conditions of the Geant4 Software License, included in the file *
// * LICENSE and available at http://cern.ch/geant4/license . These *
// * include a list of copyright holders. *
// * *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work make any representation or warranty, express or implied, *
// * regarding this software system or assume any liability for its *
// * use. Please see the license in the file LICENSE and URL above *
// * for the full disclaimer and the limitation of liability. *
// * *
// * This code implementation is the result of the scientific and *
// * technical work of the GEANT4 collaboration. *
// * By using, copying, modifying or distributing the software (or *
// * any work based on the software) you agree to acknowledge its *
// * use in resulting scientific publications, and indicate your *
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
// $Id: G4AdjointIonIonisationModel.cc,v 1.2 2009/11/20 10:31:20 ldesorgh Exp $
// GEANT4 tag $Name: geant4-09-03 $
//
#include "G4AdjointIonIonisationModel.hh"
#include "G4AdjointCSManager.hh"
#include "G4Integrator.hh"
#include "G4TrackStatus.hh"
#include "G4ParticleChange.hh"
#include "G4AdjointElectron.hh"
#include "G4AdjointProton.hh"
#include "G4AdjointInterpolator.hh"
#include "G4BetheBlochModel.hh"
#include "G4BraggIonModel.hh"
#include "G4Proton.hh"
#include "G4GenericIon.hh"
#include "G4NistManager.hh"
////////////////////////////////////////////////////////////////////////////////
//
G4AdjointIonIonisationModel::G4AdjointIonIonisationModel():
G4VEmAdjointModel("Adjoint_IonIonisation")
{
UseMatrix =true;
UseMatrixPerElement = true;
ApplyCutInRange = true;
UseOnlyOneMatrixForAllElements = true;
CS_biasing_factor =1.;
second_part_of_same_type =false;
use_only_bragg = false; // for the Ion ionisation using the parametrised table model the cross sections and the sample of secondaries is done
// as in the BraggIonModel, Therefore the use of this flag;
//The direct EM Model is taken has BetheBloch it is only used for the computation
// of the differential cross section.
//The Bragg model could be used as an alternative as it offers the same differential cross section
theBetheBlochDirectEMModel = new G4BetheBlochModel(G4GenericIon::GenericIon());
theBraggIonDirectEMModel = new G4BraggIonModel(G4GenericIon::GenericIon());
theAdjEquivOfDirectSecondPartDef=G4AdjointElectron::AdjointElectron();
theDirectPrimaryPartDef =0;
theAdjEquivOfDirectPrimPartDef =0;
/* theDirectPrimaryPartDef =fwd_ion;
theAdjEquivOfDirectPrimPartDef =adj_ion;
DefineProjectileProperty();*/
}
////////////////////////////////////////////////////////////////////////////////
//
G4AdjointIonIonisationModel::~G4AdjointIonIonisationModel()
{;}
////////////////////////////////////////////////////////////////////////////////
//
void G4AdjointIonIonisationModel::SampleSecondaries(const G4Track& aTrack,
G4bool IsScatProjToProjCase,
G4ParticleChange* fParticleChange)
{
const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
//Elastic inverse scattering
//---------------------------------------------------------
G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
return;
}
//Sample secondary energy
//-----------------------
G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, IsScatProjToProjCase);
CorrectPostStepWeight(fParticleChange,
aTrack.GetWeight(),
adjointPrimKinEnergy,
projectileKinEnergy,
IsScatProjToProjCase); //Caution !!!this weight correction should be always applied
//Kinematic:
//we consider a two body elastic scattering for the forward processes where the projectile knock on an e- at rest and gives
// him part of its energy
//----------------------------------------------------------------------------------------
G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
//Companion
//-----------
G4double companionM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
if (IsScatProjToProjCase) {
companionM0=theAdjEquivOfDirectSecondPartDef->GetPDGMass();
}
G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;
//Projectile momentum
//--------------------
G4double P_parallel = (adjointPrimP*adjointPrimP + projectileP2 - companionP2)/(2.*adjointPrimP);
G4double P_perp = std::sqrt( projectileP2 - P_parallel*P_parallel);
G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
G4double phi =G4UniformRand()*2.*3.1415926;
G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
projectileMomentum.rotateUz(dir_parallel);
if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
fParticleChange->ProposeTrackStatus(fStopAndKill);
fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
//G4cout<<"projectileMomentum "<<projectileMomentum<<G4endl;
}
else {
fParticleChange->ProposeEnergy(projectileKinEnergy);
fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
}
}
////////////////////////////////////////////////////////////////////////////////
//
G4double G4AdjointIonIonisationModel::DiffCrossSectionPerAtomPrimToSecond(
G4double kinEnergyProj,
G4double kinEnergyProd,
G4double Z,
G4double A)
{//Probably that here the Bragg Model should be also used for kinEnergyProj/nuc<2MeV
G4double dSigmadEprod=0;
G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
G4double kinEnergyProjScaled = massRatio*kinEnergyProj;
if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ //the produced particle should have a kinetic energy smaller than the projectile
G4double Tmax=kinEnergyProj;
G4double E1=kinEnergyProd;
G4double E2=kinEnergyProd*1.000001;
G4double dE=(E2-E1);
G4double sigma1,sigma2;
theDirectEMModel =theBraggIonDirectEMModel;
if (kinEnergyProjScaled >2.*MeV && !use_only_bragg) theDirectEMModel = theBetheBlochDirectEMModel; //Bethe Bloch Model
sigma1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
sigma2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
dSigmadEprod=(sigma1-sigma2)/dE;
//G4double chargeSqRatio = currentModel->GetChargeSquareRatio(theDirectPrimaryPartDef,currentMaterial,E);
if (dSigmadEprod>1.) {
G4cout<<"sigma1 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma1<<G4endl;
G4cout<<"sigma2 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma2<<G4endl;
G4cout<<"dsigma "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<dSigmadEprod<<G4endl;
}
if (theDirectEMModel == theBetheBlochDirectEMModel ){
//correction of differential cross section at high energy to correct for the suppression of particle at secondary at high
//energy used in the Bethe Bloch Model. This correction consist to multiply by g the probability function used
//to test the rejection of a secondary
//-------------------------
//Source code taken from G4BetheBlochModel::SampleSecondaries
G4double deltaKinEnergy = kinEnergyProd;
//Part of the taken code
//----------------------
// projectile formfactor - suppresion of high energy
// delta-electron production at high energy
G4double x = formfact*deltaKinEnergy;
if(x > 1.e-6) {
G4double totEnergy = kinEnergyProj + mass;
G4double etot2 = totEnergy*totEnergy;
G4double beta2 = kinEnergyProj*(kinEnergyProj + 2.0*mass)/etot2;
G4double f;
G4double f1 = 0.0;
f = 1.0 - beta2*deltaKinEnergy/Tmax;
if( 0.5 == spin ) {
f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
f += f1;
}
G4double x1 = 1.0 + x;
G4double g = 1.0/(x1*x1);
if( 0.5 == spin ) {
G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
g *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
}
if(g > 1.0) {
G4cout << "### G4BetheBlochModel in Adjoint Sim WARNING: g= " << g
<< G4endl;
g=1.;
}
//G4cout<<"g"<<g<<G4endl;
dSigmadEprod*=g;
}
}
}
return dSigmadEprod;
}
//////////////////////////////////////////////////////////////////////////////////////////////
//
void G4AdjointIonIonisationModel::SetIon(G4ParticleDefinition* adj_ion, G4ParticleDefinition* fwd_ion)
{ theDirectPrimaryPartDef =fwd_ion;
theAdjEquivOfDirectPrimPartDef =adj_ion;
DefineProjectileProperty();
}
//////////////////////////////////////////////////////////////////////////////////////////////
//
void G4AdjointIonIonisationModel::CorrectPostStepWeight(G4ParticleChange* fParticleChange, G4double old_weight,
G4double adjointPrimKinEnergy, G4double projectileKinEnergy,G4bool )
{
//It is needed because the direct cross section used to compute the differential cross section is not the one used in
// the direct model where the GenericIon stuff is considered with correction of effective charge. In the direct model the samnepl of secondaries does
// not reflect the integral cross section. The integral fwd cross section that we used to compute the differential CS
// match the sample of secondaries in the forward case despite the fact that its is not the same total CS than in the FWD case. For this reasion an extra
// weight correction is needed at the end.
G4double new_weight=old_weight;
//the correction of CS due to the problem explained above
G4double kinEnergyProjScaled = massRatio*projectileKinEnergy;
theDirectEMModel =theBraggIonDirectEMModel;
if (kinEnergyProjScaled >2.*MeV && !use_only_bragg) theDirectEMModel = theBetheBlochDirectEMModel; //Bethe Bloch Model
G4double UsedFwdCS=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,projectileKinEnergy,1,1 ,currentTcutForDirectSecond,1.e20);
G4double chargeSqRatio =1.;
if (chargeSquare>1.) chargeSqRatio = theDirectEMModel->GetChargeSquareRatio(theDirectPrimaryPartDef,currentMaterial,projectileKinEnergy);
G4double CorrectFwdCS = chargeSqRatio*theDirectEMModel->ComputeCrossSectionPerAtom(G4GenericIon::GenericIon(),kinEnergyProjScaled,1,1 ,currentTcutForDirectSecond,1.e20);
if (UsedFwdCS >0) new_weight*= CorrectFwdCS/UsedFwdCS;//May be some check is needed if UsedFwdCS ==0 probably that then we should avoid a secondary to be produced,
//additional CS crorrection needed for cross section biasing in general.
//May be wrong for ions!!! Most of the time not used!
G4double w_corr =1./CS_biasing_factor;
w_corr*=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
new_weight*=w_corr;
new_weight*=projectileKinEnergy/adjointPrimKinEnergy;
fParticleChange->SetParentWeightByProcess(false);
fParticleChange->SetSecondaryWeightByProcess(false);
fParticleChange->ProposeParentWeight(new_weight);
}
//////////////////////////////////////////////////////////////////////////////////////////////
//
void G4AdjointIonIonisationModel::DefineProjectileProperty()
{
//Slightly modified code taken from G4BetheBlochModel::SetParticle
//------------------------------------------------
G4String pname = theDirectPrimaryPartDef->GetParticleName();
if (theDirectPrimaryPartDef->GetParticleType() == "nucleus" &&
pname != "deuteron" && pname != "triton") {
isIon = true;
}
mass = theDirectPrimaryPartDef->GetPDGMass();
massRatio= G4GenericIon::GenericIon()->GetPDGMass()/mass;
spin = theDirectPrimaryPartDef->GetPDGSpin();
G4double q = theDirectPrimaryPartDef->GetPDGCharge()/eplus;
chargeSquare = q*q;
ratio = electron_mass_c2/mass;
ratio2 = ratio*ratio;
one_plus_ratio_2=(1+ratio)*(1+ratio);
one_minus_ratio_2=(1-ratio)*(1-ratio);
G4double magmom = theDirectPrimaryPartDef->GetPDGMagneticMoment()
*mass/(0.5*eplus*hbar_Planck*c_squared);
magMoment2 = magmom*magmom - 1.0;
formfact = 0.0;
if(theDirectPrimaryPartDef->GetLeptonNumber() == 0) {
G4double x = 0.8426*GeV;
if(spin == 0.0 && mass < GeV) {x = 0.736*GeV;}
else if(mass > GeV) {
x /= G4NistManager::Instance()->GetZ13(mass/proton_mass_c2);
// tlimit = 51.2*GeV*A13[iz]*A13[iz];
}
formfact = 2.0*electron_mass_c2/(x*x);
tlimit = 2.0/formfact;
}
}
//////////////////////////////////////////////////////////////////////////////
//
G4double G4AdjointIonIonisationModel::GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
{
G4double Tmax=PrimAdjEnergy*one_plus_ratio_2/(one_minus_ratio_2-2.*ratio*PrimAdjEnergy/mass);
return Tmax;
}
//////////////////////////////////////////////////////////////////////////////
//
G4double G4AdjointIonIonisationModel::GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy,G4double Tcut)
{ return PrimAdjEnergy+Tcut;
}
//////////////////////////////////////////////////////////////////////////////
//
G4double G4AdjointIonIonisationModel::GetSecondAdjEnergyMaxForProdToProjCase(G4double )
{ return HighEnergyLimit;
}
//////////////////////////////////////////////////////////////////////////////
//
G4double G4AdjointIonIonisationModel::GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
{ G4double Tmin= (2*PrimAdjEnergy-4*mass + std::sqrt(4.*PrimAdjEnergy*PrimAdjEnergy +16.*mass*mass + 8.*PrimAdjEnergy*mass*(1/ratio +ratio)))/4.;
return Tmin;
}
source/processes/electromagnetic/adjoint/src/G4AdjointPhotoElectricModel.cc
+72 -38
View File
@@ -23,6 +23,9 @@
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
// $Id: G4AdjointPhotoElectricModel.cc,v 1.5 2009/12/16 17:50:05 gunter Exp $
// GEANT4 tag $Name: geant4-09-03 $
//
#include "G4AdjointPhotoElectricModel.hh"
#include "G4AdjointCSManager.hh"
@@ -34,14 +37,22 @@
#include "G4Gamma.hh"
#include "G4AdjointGamma.hh"
////////////////////////////////////////////////////////////////////////////////
//
G4AdjointPhotoElectricModel::G4AdjointPhotoElectricModel():
G4VEmAdjointModel("AdjointPEEffect")
{ SetUseMatrix(false);
SetApplyCutInRange(false);
current_eEnergy =0.;
totAdjointCS=0.;
theAdjEquivOfDirectPrimPartDef =G4AdjointGamma::AdjointGamma();
theAdjEquivOfDirectSecondPartDef=G4AdjointElectron::AdjointElectron();
theDirectPrimaryPartDef=G4Gamma::Gamma();
second_part_of_same_type=false;
theDirectPEEffectModel = new G4PEEffectModel();
}
////////////////////////////////////////////////////////////////////////////////
//
@@ -56,47 +67,32 @@ void G4AdjointPhotoElectricModel::SampleSecondaries(const G4Track& aTrack,
{ if (IsScatProjToProjCase) return ;
//Compute the totAdjointCS vectors if not already done for the current couple and electron energy
//-----------------------------------------------------------------------------------------------
const G4MaterialCutsCouple* aCouple = aTrack.GetMaterialCutsCouple();
const G4DynamicParticle* aDynPart = aTrack.GetDynamicParticle() ;
G4double electronEnergy = aDynPart->GetKineticEnergy();
G4ThreeVector electronDirection= aDynPart->GetMomentumDirection() ;
totAdjointCS = AdjointCrossSection(aCouple, electronEnergy,IsScatProjToProjCase);
pre_step_AdjointCS = totAdjointCS; //The last computed CS was at pre step point
G4double adjCS;
adjCS = AdjointCrossSection(aCouple, electronEnergy,IsScatProjToProjCase);
post_step_AdjointCS = totAdjointCS;
//Sample gamma energy
//-------------
/////////////////////////////////////////////////////////////////////////////////
// Module: G4ContinuousGainOfEnergy.hh
// Author: L. Desorgher
// Date: 1 September 2007
// Organisation: SpaceIT GmbH
// Customer: ESA/ESTEC
/////////////////////////////////////////////////////////////////////////////////
//
// CHANGE HISTORY
// --------------
// ChangeHistory:
// 1 September 2007 creation by L. Desorgher
//
//-------------------------------------------------------------
// Documentation:
// Modell for the adjoint compton scattering
//
//Sample element
//-------------
const G4ElementVector* theElementVector = currentMaterial->GetElementVector();
const G4double* theAtomNumDensityVector = currentMaterial->GetVecNbOfAtomsPerVolume();
size_t nelm = currentMaterial->GetNumberOfElements();
G4double rand_CS= totAdjointCS*G4UniformRand();
G4double rand_CS= G4UniformRand()*xsec[nelm-1];
for (index_element=0; index_element<nelm-1; index_element++){
if (rand_CS<xsec[index_element]) break;
}
//Sample shell and binding energy
//-------------
rand_CS= totAdjointCS*G4UniformRand()/theAtomNumDensityVector[index_element];
G4int nShells = (*theElementVector)[index_element]->GetNbOfAtomicShells();
rand_CS= shell_prob[index_element][nShells-1]*G4UniformRand();
G4int i = 0;
for (i=0; i<nShells-1; i++){
if (rand_CS<shell_prob[index_element][i]) break;
@@ -140,7 +136,7 @@ void G4AdjointPhotoElectricModel::SampleSecondaries(const G4Track& aTrack,
//Weight correction
//-----------------------
CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(), electronEnergy,gammaEnergy);
CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(), electronEnergy,gammaEnergy,IsScatProjToProjCase);
@@ -148,13 +144,36 @@ void G4AdjointPhotoElectricModel::SampleSecondaries(const G4Track& aTrack,
//--------------------------------------------
G4DynamicParticle* anAdjointGamma = new G4DynamicParticle (
G4AdjointGamma::AdjointGamma(),adjoint_gammaDirection, gammaEnergy);
fParticleChange->ProposeTrackStatus(fStopAndKill);
fParticleChange->AddSecondary(anAdjointGamma);
fParticleChange->AddSecondary(anAdjointGamma);
}
////////////////////////////////////////////////////////////////////////////////
//
void G4AdjointPhotoElectricModel::CorrectPostStepWeight(G4ParticleChange* fParticleChange,
G4double old_weight,
G4double adjointPrimKinEnergy,
G4double projectileKinEnergy ,
G4bool )
{
G4double new_weight=old_weight;
G4double w_corr =G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection()/factorCSBiasing;
w_corr*=post_step_AdjointCS/pre_step_AdjointCS;
new_weight*=w_corr;
new_weight*=projectileKinEnergy/adjointPrimKinEnergy;
fParticleChange->SetParentWeightByProcess(false);
fParticleChange->SetSecondaryWeightByProcess(false);
fParticleChange->ProposeParentWeight(new_weight);
}
////////////////////////////////////////////////////////////////////////////////
@@ -163,21 +182,33 @@ void G4AdjointPhotoElectricModel::SampleSecondaries(const G4Track& aTrack,
G4double G4AdjointPhotoElectricModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
G4double electronEnergy,
G4bool IsScatProjToProjCase)
{ if (IsScatProjToProjCase) return 0.;
{
if (IsScatProjToProjCase) return 0.;
if (aCouple !=currentCouple || current_eEnergy !=electronEnergy) {
totAdjointCS = 0.;
DefineCurrentMaterialAndElectronEnergy(aCouple, electronEnergy);
const G4ElementVector* theElementVector = currentMaterial->GetElementVector();
const G4double* theAtomNumDensityVector = currentMaterial->GetVecNbOfAtomsPerVolume();
const double* theAtomNumDensityVector = currentMaterial->GetVecNbOfAtomsPerVolume();
size_t nelm = currentMaterial->GetNumberOfElements();
for (index_element=0;index_element<nelm;index_element++){
totAdjointCS +=AdjointCrossSectionPerAtom((*theElementVector)[index_element],electronEnergy)*theAtomNumDensityVector[index_element];
xsec[index_element] = totAdjointCS;
}
}
return totAdjointCS;
totBiasedAdjointCS=std::min(totAdjointCS,0.01);
// totBiasedAdjointCS=totAdjointCS;
factorCSBiasing = totBiasedAdjointCS/totAdjointCS;
lastCS=totBiasedAdjointCS;
}
return totBiasedAdjointCS;
}
////////////////////////////////////////////////////////////////////////////////
@@ -187,25 +218,28 @@ G4double G4AdjointPhotoElectricModel::AdjointCrossSectionPerAtom(const G4Element
{
G4int nShells = anElement->GetNbOfAtomicShells();
G4double Z= anElement->GetZ();
G4double N= anElement->GetN();
G4int i = 0;
G4double B0=anElement->GetAtomicShell(0);
G4double gammaEnergy = electronEnergy+B0;
G4double adjointCS = theDirectPEEffectModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(),gammaEnergy,Z,N,0.,0.)*electronEnergy/gammaEnergy;
G4double CS= theDirectPEEffectModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(),gammaEnergy,Z,0.,0.,0.);
G4double adjointCS =0.;
if (CS >0) adjointCS += CS/gammaEnergy;
shell_prob[index_element][0] = adjointCS;
for (i=1;i<nShells;i++){
//G4cout<<i<<std::endl;
//G4cout<<i<<G4endl;
G4double Bi_= anElement->GetAtomicShell(i-1);
G4double Bi = anElement->GetAtomicShell(i);
//G4cout<<Bi_<<'\t'<<Bi<<std::endl;
//G4cout<<Bi_<<'\t'<<Bi<<G4endl;
if (electronEnergy <Bi_-Bi) {
gammaEnergy = electronEnergy+Bi;
adjointCS +=theDirectPEEffectModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(),gammaEnergy,anElement->GetZ(),N,0.,0.)*electronEnergy/gammaEnergy;
CS=theDirectPEEffectModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(),gammaEnergy,Z,0.,0.,0.);
if (CS>0) adjointCS +=CS/gammaEnergy;
}
shell_prob[index_element][i] = adjointCS;
}
adjointCS*=electronEnergy;
return adjointCS;
}
source/processes/electromagnetic/adjoint/src/G4AdjointProcessEquivalentToDirectProcess.cc
+321
View File
@@ -0,0 +1,321 @@
//
// ********************************************************************
// * License and Disclaimer *
// * *
// * The Geant4 software is copyright of the Copyright Holders of *
// * the Geant4 Collaboration. It is provided under the terms and *
// * conditions of the Geant4 Software License, included in the file *
// * LICENSE and available at http://cern.ch/geant4/license . These *
// * include a list of copyright holders. *
// * *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work make any representation or warranty, express or implied, *
// * regarding this software system or assume any liability for its *
// * use. Please see the license in the file LICENSE and URL above *
// * for the full disclaimer and the limitation of liability. *
// * *
// * This code implementation is the result of the scientific and *
// * technical work of the GEANT4 collaboration. *
// * By using, copying, modifying or distributing the software (or *
// * any work based on the software) you agree to acknowledge its *
// * use in resulting scientific publications, and indicate your *
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
//
// $Id: G4AdjointProcessEquivalentToDirectProcess.cc,v 1.1 2009/11/11 00:31:19 ldesorgh Exp $
// GEANT4 tag $Name: geant4-09-03 $
//
//
// ------------------------------------------------------------
// GEANT 4 class implementation file
//
// Class Description
//
// This class is for adjoint process equivalent to direct process
// ------------------------------------------------------------
// Created by L.Desorgher 25 Sept. 2009 Inspired from G4WrapperProcess
// ------------------------------------------------------------
#include "G4AdjointProcessEquivalentToDirectProcess.hh"
#include "G4DynamicParticle.hh"
G4AdjointProcessEquivalentToDirectProcess::G4AdjointProcessEquivalentToDirectProcess(const G4String& aName,
G4VProcess* aProcess,
G4ParticleDefinition* fwd_particle_def)
:G4VProcess(aName)
{
theDirectProcess =aProcess;
theProcessType = theDirectProcess->GetProcessType();
theFwdParticleDef = fwd_particle_def;
}
G4AdjointProcessEquivalentToDirectProcess::~G4AdjointProcessEquivalentToDirectProcess()
{
if (theDirectProcess!=0) delete theDirectProcess;
}
void G4AdjointProcessEquivalentToDirectProcess::ResetNumberOfInteractionLengthLeft()
{
theDirectProcess->ResetNumberOfInteractionLengthLeft();
}
G4double G4AdjointProcessEquivalentToDirectProcess::
AlongStepGetPhysicalInteractionLength( const G4Track& track,
G4double previousStepSize,
G4double currentMinimumStep,
G4double& proposedSafety,
G4GPILSelection* selection )
{
//Change the particle definition to the direct one
//------------------------------------------------
G4DynamicParticle* theDynPart = const_cast<G4DynamicParticle*> (track.GetDynamicParticle());
G4ParticleDefinition* adjPartDef = theDynPart->GetDefinition();
G4DecayProducts* decayProducts = const_cast<G4DecayProducts*> (theDynPart->GetPreAssignedDecayProducts());
theDynPart->SetPreAssignedDecayProducts((G4DecayProducts*)(0));
theDynPart->SetDefinition(theFwdParticleDef);
//Call the direct process
//----------------------
G4double GPIL = theDirectProcess->
AlongStepGetPhysicalInteractionLength( track,
previousStepSize,
currentMinimumStep,
proposedSafety,
selection );
//Restore the adjoint particle definition to the direct one
//------------------------------------------------
theDynPart->SetDefinition(adjPartDef);
theDynPart->SetPreAssignedDecayProducts(decayProducts);
return GPIL;
}
G4double G4AdjointProcessEquivalentToDirectProcess::
AtRestGetPhysicalInteractionLength( const G4Track& track,
G4ForceCondition* condition )
{ //Change the particle definition to the direct one
//------------------------------------------------
G4DynamicParticle* theDynPart = const_cast<G4DynamicParticle*> (track.GetDynamicParticle());
G4ParticleDefinition* adjPartDef = theDynPart->GetDefinition();
G4DecayProducts* decayProducts = const_cast<G4DecayProducts*> (theDynPart->GetPreAssignedDecayProducts());
theDynPart->SetPreAssignedDecayProducts((G4DecayProducts*)(0));
theDynPart->SetDefinition(theFwdParticleDef);
//Call the direct process
//----------------------
G4double GPIL = theDirectProcess->AtRestGetPhysicalInteractionLength( track, condition );
//Restore the adjoint particle definition to the direct one
//------------------------------------------------
theDynPart->SetDefinition(adjPartDef);
theDynPart->SetPreAssignedDecayProducts(decayProducts);
return GPIL;
}
G4double G4AdjointProcessEquivalentToDirectProcess::
PostStepGetPhysicalInteractionLength( const G4Track& track,
G4double previousStepSize,
G4ForceCondition* condition )
{
//Change the particle definition to the direct one
//------------------------------------------------
G4DynamicParticle* theDynPart = const_cast<G4DynamicParticle*> (track.GetDynamicParticle());
G4ParticleDefinition* adjPartDef = theDynPart->GetDefinition();
G4DecayProducts* decayProducts = const_cast<G4DecayProducts*> (theDynPart->GetPreAssignedDecayProducts());
theDynPart->SetPreAssignedDecayProducts((G4DecayProducts*)(0));
theDynPart->SetDefinition(theFwdParticleDef);
//Call the direct process
//----------------------
G4double GPIL = theDirectProcess->PostStepGetPhysicalInteractionLength( track,
previousStepSize,
condition );
//Restore the adjoint particle definition to the direct one
//------------------------------------------------
theDynPart->SetDefinition(adjPartDef);
theDynPart->SetPreAssignedDecayProducts(decayProducts);
return GPIL;
}
/*
void G4AdjointProcessEquivalentToDirectProcess::SetProcessManager(const G4ProcessManager* procMan)
{
theDirectProcess->SetProcessManager(procMan);
}
const G4ProcessManager* G4AdjointProcessEquivalentToDirectProcess::GetProcessManager()
{
return theDirectProcess->GetProcessManager();
}
*/
G4VParticleChange* G4AdjointProcessEquivalentToDirectProcess::PostStepDoIt( const G4Track& track,
const G4Step& stepData )
{
//Change the particle definition to the direct one
//------------------------------------------------
G4DynamicParticle* theDynPart = const_cast<G4DynamicParticle*> (track.GetDynamicParticle());
G4ParticleDefinition* adjPartDef = theDynPart->GetDefinition();
G4DecayProducts* decayProducts = const_cast<G4DecayProducts*> (theDynPart->GetPreAssignedDecayProducts());
theDynPart->SetPreAssignedDecayProducts((G4DecayProducts*)(0));
theDynPart->SetDefinition(theFwdParticleDef);
//Call the direct process
//----------------------
G4VParticleChange* partChange = theDirectProcess->PostStepDoIt( track, stepData );
//Restore the adjoint particle definition to the direct one
//------------------------------------------------
theDynPart->SetDefinition(adjPartDef);
theDynPart->SetPreAssignedDecayProducts(decayProducts);
return partChange;
}
G4VParticleChange* G4AdjointProcessEquivalentToDirectProcess::AlongStepDoIt( const G4Track& track,
const G4Step& stepData )
{
//Change the particle definition to the direct one
//------------------------------------------------
G4DynamicParticle* theDynPart = const_cast<G4DynamicParticle*> (track.GetDynamicParticle());
G4ParticleDefinition* adjPartDef = theDynPart->GetDefinition();
G4DecayProducts* decayProducts = const_cast<G4DecayProducts*> (theDynPart->GetPreAssignedDecayProducts());
theDynPart->SetPreAssignedDecayProducts((G4DecayProducts*)(0));
theDynPart->SetDefinition(theFwdParticleDef);
//Call the direct process
//----------------------
G4VParticleChange* partChange =theDirectProcess->AlongStepDoIt( track, stepData );
//Restore the adjoint particle definition to the direct one
//------------------------------------------------
theDynPart->SetDefinition(adjPartDef);
theDynPart->SetPreAssignedDecayProducts(decayProducts);
return partChange;
}
G4VParticleChange* G4AdjointProcessEquivalentToDirectProcess::AtRestDoIt( const G4Track& track,
const G4Step& stepData )
{
//Change the particle definition to the direct one
//------------------------------------------------
G4DynamicParticle* theDynPart = const_cast<G4DynamicParticle*> (track.GetDynamicParticle());
G4ParticleDefinition* adjPartDef = theDynPart->GetDefinition();
G4DecayProducts* decayProducts = const_cast<G4DecayProducts*> (theDynPart->GetPreAssignedDecayProducts());
theDynPart->SetPreAssignedDecayProducts((G4DecayProducts*)(0));
theDynPart->SetDefinition(theFwdParticleDef);
//Call the direct process
//----------------------
G4VParticleChange* partChange =theDirectProcess->AtRestDoIt( track, stepData );
//Restore the adjoint particle definition to the direct one
//------------------------------------------------
theDynPart->SetDefinition(adjPartDef);
theDynPart->SetPreAssignedDecayProducts(decayProducts);
return partChange;
}
G4bool G4AdjointProcessEquivalentToDirectProcess::IsApplicable(const G4ParticleDefinition&)
{
return theDirectProcess->IsApplicable(*theFwdParticleDef);
}
void G4AdjointProcessEquivalentToDirectProcess::BuildPhysicsTable(const G4ParticleDefinition& )
{
return theDirectProcess->BuildPhysicsTable(*theFwdParticleDef);
}
void G4AdjointProcessEquivalentToDirectProcess::PreparePhysicsTable(const G4ParticleDefinition& )
{
return theDirectProcess->PreparePhysicsTable(*theFwdParticleDef);
}
G4bool G4AdjointProcessEquivalentToDirectProcess::
StorePhysicsTable(const G4ParticleDefinition* ,
const G4String& directory,
G4bool ascii)
{
return theDirectProcess->StorePhysicsTable(theFwdParticleDef, directory, ascii);
}
G4bool G4AdjointProcessEquivalentToDirectProcess::
RetrievePhysicsTable( const G4ParticleDefinition* ,
const G4String& directory,
G4bool ascii)
{
return theDirectProcess->RetrievePhysicsTable(theFwdParticleDef, directory, ascii);
}
void G4AdjointProcessEquivalentToDirectProcess::StartTracking(G4Track* track)
{
//Change the particle definition to the direct one
//------------------------------------------------
G4DynamicParticle* theDynPart = const_cast<G4DynamicParticle*> (track->GetDynamicParticle());
G4ParticleDefinition* adjPartDef = theDynPart->GetDefinition();
G4DecayProducts* decayProducts = const_cast<G4DecayProducts*> (theDynPart->GetPreAssignedDecayProducts());
theDynPart->SetPreAssignedDecayProducts((G4DecayProducts*)(0));
theDynPart->SetDefinition(theFwdParticleDef);
theDirectProcess->StartTracking(track);
//Restore the adjoint particle definition to the direct one
//------------------------------------------------
theDynPart->SetDefinition(adjPartDef);
theDynPart->SetPreAssignedDecayProducts(decayProducts);
return;
}
void G4AdjointProcessEquivalentToDirectProcess::EndTracking()
{
theDirectProcess->EndTracking();
}
source/processes/electromagnetic/adjoint/src/G4AdjointeIonisationModel.cc
+216
View File
@@ -0,0 +1,216 @@
//
// ********************************************************************
// * License and Disclaimer *
// * *
// * The Geant4 software is copyright of the Copyright Holders of *
// * the Geant4 Collaboration. It is provided under the terms and *
// * conditions of the Geant4 Software License, included in the file *
// * LICENSE and available at http://cern.ch/geant4/license . These *
// * include a list of copyright holders. *
// * *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work make any representation or warranty, express or implied, *
// * regarding this software system or assume any liability for its *
// * use. Please see the license in the file LICENSE and URL above *
// * for the full disclaimer and the limitation of liability. *
// * *
// * This code implementation is the result of the scientific and *
// * technical work of the GEANT4 collaboration. *
// * By using, copying, modifying or distributing the software (or *
// * any work based on the software) you agree to acknowledge its *
// * use in resulting scientific publications, and indicate your *
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
// $Id: G4AdjointeIonisationModel.cc,v 1.2 2009/11/20 10:31:20 ldesorgh Exp $
// GEANT4 tag $Name: geant4-09-03 $
//
#include "G4AdjointeIonisationModel.hh"
#include "G4AdjointCSManager.hh"
#include "G4Integrator.hh"
#include "G4TrackStatus.hh"
#include "G4ParticleChange.hh"
#include "G4AdjointElectron.hh"
#include "G4Gamma.hh"
#include "G4AdjointGamma.hh"
////////////////////////////////////////////////////////////////////////////////
//
G4AdjointeIonisationModel::G4AdjointeIonisationModel():
G4VEmAdjointModel("Inv_eIon_model")
{
UseMatrix =true;
UseMatrixPerElement = true;
ApplyCutInRange = true;
UseOnlyOneMatrixForAllElements = true;
CS_biasing_factor =1.;
WithRapidSampling = false;
theAdjEquivOfDirectPrimPartDef =G4AdjointElectron::AdjointElectron();
theAdjEquivOfDirectSecondPartDef=G4AdjointElectron::AdjointElectron();
theDirectPrimaryPartDef=G4Electron::Electron();
second_part_of_same_type=true;
}
////////////////////////////////////////////////////////////////////////////////
//
G4AdjointeIonisationModel::~G4AdjointeIonisationModel()
{;}
////////////////////////////////////////////////////////////////////////////////
//
void G4AdjointeIonisationModel::SampleSecondaries(const G4Track& aTrack,
G4bool IsScatProjToProjCase,
G4ParticleChange* fParticleChange)
{
const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
//Elastic inverse scattering
//---------------------------------------------------------
G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
return;
}
//Sample secondary energy
//-----------------------
G4double projectileKinEnergy;
if (!WithRapidSampling ) { //used by default
projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, IsScatProjToProjCase);
CorrectPostStepWeight(fParticleChange,
aTrack.GetWeight(),
adjointPrimKinEnergy,
projectileKinEnergy,
IsScatProjToProjCase); //Caution !!!this weight correction should be always applied
}
else { //only for test at the moment
G4double Emin,Emax;
if (IsScatProjToProjCase) {
Emin=GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond);
Emax=GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
}
else {
Emin=GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);
Emax=GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
}
projectileKinEnergy = Emin*std::pow(Emax/Emin,G4UniformRand());
lastCS=lastAdjointCSForScatProjToProjCase;
if ( !IsScatProjToProjCase) lastCS=lastAdjointCSForProdToProjCase;
G4double new_weight=aTrack.GetWeight();
G4double used_diffCS=lastCS*std::log(Emax/Emin)/projectileKinEnergy;
G4double needed_diffCS=adjointPrimKinEnergy/projectileKinEnergy;
if (!IsScatProjToProjCase) needed_diffCS *=DiffCrossSectionPerVolumePrimToSecond(currentMaterial,projectileKinEnergy,adjointPrimKinEnergy);
else needed_diffCS *=DiffCrossSectionPerVolumePrimToScatPrim(currentMaterial,projectileKinEnergy,adjointPrimKinEnergy);
new_weight*=needed_diffCS/used_diffCS;
fParticleChange->SetParentWeightByProcess(false);
fParticleChange->SetSecondaryWeightByProcess(false);
fParticleChange->ProposeParentWeight(new_weight);
}
//Kinematic:
//we consider a two body elastic scattering for the forward processes where the projectile knock on an e- at rest and gives
// him part of its energy
//----------------------------------------------------------------------------------------
G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
//Companion
//-----------
G4double companionM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
if (IsScatProjToProjCase) {
companionM0=theAdjEquivOfDirectSecondPartDef->GetPDGMass();
}
G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;
//Projectile momentum
//--------------------
G4double P_parallel = (adjointPrimP*adjointPrimP + projectileP2 - companionP2)/(2.*adjointPrimP);
G4double P_perp = std::sqrt( projectileP2 - P_parallel*P_parallel);
G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
G4double phi =G4UniformRand()*2.*3.1415926;
G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
projectileMomentum.rotateUz(dir_parallel);
if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
fParticleChange->ProposeTrackStatus(fStopAndKill);
fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
//G4cout<<"projectileMomentum "<<projectileMomentum<<G4endl;
}
else {
fParticleChange->ProposeEnergy(projectileKinEnergy);
fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
}
}
////////////////////////////////////////////////////////////////////////////////
//
//The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
G4double G4AdjointeIonisationModel::DiffCrossSectionPerAtomPrimToSecond(
G4double kinEnergyProj,
G4double kinEnergyProd,
G4double Z,
G4double )
{
G4double dSigmadEprod=0;
G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ //the produced particle should have a kinetic energy smaller than the projectile
dSigmadEprod=Z*DiffCrossSectionMoller(kinEnergyProj,kinEnergyProd);
}
return dSigmadEprod;
}
//////////////////////////////////////////////////////////////////////////////
//
G4double G4AdjointeIonisationModel::DiffCrossSectionMoller(G4double kinEnergyProj,G4double kinEnergyProd){
G4double electron_mass_c2=0.51099906*MeV;
G4double energy = kinEnergyProj + electron_mass_c2;
G4double x = kinEnergyProd/kinEnergyProj;
G4double gam = energy/electron_mass_c2;
G4double gamma2 = gam*gam;
G4double beta2 = 1.0 - 1.0/gamma2;
G4double g = (2.0*gam - 1.0)/gamma2;
G4double y = 1.0 - x;
G4double fac=twopi_mc2_rcl2/electron_mass_c2;
G4double dCS = fac*( 1.-g + ((1.0 - g*x)/(x*x)) + ((1.0 - g*y)/(y*y)))/(beta2*(gam-1));
return dCS/kinEnergyProj;
}
source/processes/electromagnetic/adjoint/src/G4AdjointhIonisationModel.cc
+502
View File
@@ -0,0 +1,502 @@
//
// ********************************************************************
// * License and Disclaimer *
// * *
// * The Geant4 software is copyright of the Copyright Holders of *
// * the Geant4 Collaboration. It is provided under the terms and *
// * conditions of the Geant4 Software License, included in the file *
// * LICENSE and available at http://cern.ch/geant4/license . These *
// * include a list of copyright holders. *
// * *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work make any representation or warranty, express or implied, *
// * regarding this software system or assume any liability for its *
// * use. Please see the license in the file LICENSE and URL above *
// * for the full disclaimer and the limitation of liability. *
// * *
// * This code implementation is the result of the scientific and *
// * technical work of the GEANT4 collaboration. *
// * By using, copying, modifying or distributing the software (or *
// * any work based on the software) you agree to acknowledge its *
// * use in resulting scientific publications, and indicate your *
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
// $Id: G4AdjointhIonisationModel.cc,v 1.3 2009/12/16 17:50:07 gunter Exp $
// GEANT4 tag $Name: geant4-09-03 $
//
#include "G4AdjointhIonisationModel.hh"
#include "G4AdjointCSManager.hh"
#include "G4Integrator.hh"
#include "G4TrackStatus.hh"
#include "G4ParticleChange.hh"
#include "G4AdjointElectron.hh"
#include "G4AdjointProton.hh"
#include "G4AdjointInterpolator.hh"
#include "G4BetheBlochModel.hh"
#include "G4BraggModel.hh"
#include "G4Proton.hh"
#include "G4NistManager.hh"
////////////////////////////////////////////////////////////////////////////////
//
G4AdjointhIonisationModel::G4AdjointhIonisationModel(G4ParticleDefinition* projectileDefinition):
G4VEmAdjointModel("Adjoint_hIonisation")
{
UseMatrix =true;
UseMatrixPerElement = true;
ApplyCutInRange = true;
UseOnlyOneMatrixForAllElements = true;
CS_biasing_factor =1.;
second_part_of_same_type =false;
//The direct EM Modfel is taken has BetheBloch it is only used for the computation
// of the differential cross section.
//The Bragg model could be used as an alternative as it offers the same differential cross section
theDirectEMModel = new G4BetheBlochModel(projectileDefinition);
theBraggDirectEMModel = new G4BraggModel(projectileDefinition);
theAdjEquivOfDirectSecondPartDef=G4AdjointElectron::AdjointElectron();
theDirectPrimaryPartDef = projectileDefinition;
if (projectileDefinition == G4Proton::Proton()) {
theAdjEquivOfDirectPrimPartDef = G4AdjointProton::AdjointProton();
}
DefineProjectileProperty();
}
////////////////////////////////////////////////////////////////////////////////
//
G4AdjointhIonisationModel::~G4AdjointhIonisationModel()
{;}
////////////////////////////////////////////////////////////////////////////////
//
void G4AdjointhIonisationModel::SampleSecondaries(const G4Track& aTrack,
G4bool IsScatProjToProjCase,
G4ParticleChange* fParticleChange)
{
if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange);
const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
//Elastic inverse scattering
//---------------------------------------------------------
G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
return;
}
//Sample secondary energy
//-----------------------
G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, IsScatProjToProjCase);
CorrectPostStepWeight(fParticleChange,
aTrack.GetWeight(),
adjointPrimKinEnergy,
projectileKinEnergy,
IsScatProjToProjCase); //Caution !!!this weight correction should be always applied
//Kinematic:
//we consider a two body elastic scattering for the forward processes where the projectile knock on an e- at rest and gives
// him part of its energy
//----------------------------------------------------------------------------------------
G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
//Companion
//-----------
G4double companionM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
if (IsScatProjToProjCase) {
companionM0=theAdjEquivOfDirectSecondPartDef->GetPDGMass();
}
G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;
//Projectile momentum
//--------------------
G4double P_parallel = (adjointPrimP*adjointPrimP + projectileP2 - companionP2)/(2.*adjointPrimP);
G4double P_perp = std::sqrt( projectileP2 - P_parallel*P_parallel);
G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
G4double phi =G4UniformRand()*2.*3.1415926;
G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
projectileMomentum.rotateUz(dir_parallel);
if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
fParticleChange->ProposeTrackStatus(fStopAndKill);
fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
//G4cout<<"projectileMomentum "<<projectileMomentum<<G4endl;
}
else {
fParticleChange->ProposeEnergy(projectileKinEnergy);
fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
}
}
////////////////////////////////////////////////////////////////////////////////
//
void G4AdjointhIonisationModel::RapidSampleSecondaries(const G4Track& aTrack,
G4bool IsScatProjToProjCase,
G4ParticleChange* fParticleChange)
{
const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
return;
}
G4double projectileKinEnergy =0.;
G4double eEnergy=0.;
G4double newCS=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2*mass;
if (!IsScatProjToProjCase){//1/E^2 distribution
eEnergy=adjointPrimKinEnergy;
G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
G4double Emin= GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);
if (Emin>=Emax) return;
G4double a=1./Emax;
G4double b=1./Emin;
newCS=newCS*(b-a)/eEnergy;
projectileKinEnergy =1./(b- (b-a)*G4UniformRand());
}
else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
G4double Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond);
if (Emin>=Emax) return;
G4double diff1=Emin-adjointPrimKinEnergy;
G4double diff2=Emax-adjointPrimKinEnergy;
G4double t1=adjointPrimKinEnergy*(1./diff1-1./diff2);
G4double t2=adjointPrimKinEnergy*(1./Emin-1./Emax);
G4double f31=diff1/Emin;
G4double f32=diff2/Emax/f31;
G4double t3=2.*std::log(f32);
G4double sum_t=t1+t2+t3;
newCS=newCS*sum_t/adjointPrimKinEnergy/adjointPrimKinEnergy;
G4double t=G4UniformRand()*sum_t;
if (t <=t1 ){
G4double q= G4UniformRand()*t1/adjointPrimKinEnergy ;
projectileKinEnergy =adjointPrimKinEnergy +1./(1./diff1-q);
}
else if (t <=t2 ) {
G4double q= G4UniformRand()*t2/adjointPrimKinEnergy;
projectileKinEnergy =1./(1./Emin-q);
}
else {
projectileKinEnergy=adjointPrimKinEnergy/(1.-f31*std::pow(f32,G4UniformRand()));
}
eEnergy=projectileKinEnergy-adjointPrimKinEnergy;
}
G4double diffCS_perAtom_Used=twopi_mc2_rcl2*mass*adjointPrimKinEnergy/projectileKinEnergy/projectileKinEnergy/eEnergy/eEnergy;
//Weight correction
//-----------------------
//First w_corr is set to the ratio between adjoint total CS and fwd total CS
G4double w_corr=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
//G4cout<<w_corr<<G4endl;
w_corr*=newCS/lastCS;
//G4cout<<w_corr<<G4endl;
//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
//Here we consider the true diffCS as the one obtained by the numerical differentiation over Tcut of the direct CS
G4double diffCS = DiffCrossSectionPerAtomPrimToSecond(projectileKinEnergy, eEnergy,1,1);
w_corr*=diffCS/diffCS_perAtom_Used;
//G4cout<<w_corr<<G4endl;
G4double new_weight = aTrack.GetWeight()*w_corr;
fParticleChange->SetParentWeightByProcess(false);
fParticleChange->SetSecondaryWeightByProcess(false);
fParticleChange->ProposeParentWeight(new_weight);
//Kinematic:
//we consider a two body elastic scattering for the forward processes where the projectile knock on an e- at rest and gives
// him part of its energy
//----------------------------------------------------------------------------------------
G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
//Companion
//-----------
G4double companionM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
if (IsScatProjToProjCase) {
companionM0=theAdjEquivOfDirectSecondPartDef->GetPDGMass();
}
G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;
//Projectile momentum
//--------------------
G4double P_parallel = (adjointPrimP*adjointPrimP + projectileP2 - companionP2)/(2.*adjointPrimP);
G4double P_perp = std::sqrt( projectileP2 - P_parallel*P_parallel);
G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
G4double phi =G4UniformRand()*2.*3.1415926;
G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
projectileMomentum.rotateUz(dir_parallel);
if (!IsScatProjToProjCase ){ //kill the primary and add a secondary
fParticleChange->ProposeTrackStatus(fStopAndKill);
fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
//G4cout<<"projectileMomentum "<<projectileMomentum<<G4endl;
}
else {
fParticleChange->ProposeEnergy(projectileKinEnergy);
fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
}
}
////////////////////////////////////////////////////////////////////////////////
//
G4double G4AdjointhIonisationModel::DiffCrossSectionPerAtomPrimToSecond(
G4double kinEnergyProj,
G4double kinEnergyProd,
G4double Z,
G4double A)
{//Probably that here the Bragg Model should be also used for kinEnergyProj/nuc<2MeV
G4double dSigmadEprod=0;
G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ //the produced particle should have a kinetic energy smaller than the projectile
G4double Tmax=kinEnergyProj;
G4double E1=kinEnergyProd;
G4double E2=kinEnergyProd*1.000001;
G4double dE=(E2-E1);
G4double sigma1,sigma2;
if (kinEnergyProj >2.*MeV){
sigma1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
sigma2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
}
else {
sigma1=theBraggDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
sigma2=theBraggDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
}
dSigmadEprod=(sigma1-sigma2)/dE;
if (dSigmadEprod>1.) {
G4cout<<"sigma1 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma1<<G4endl;
G4cout<<"sigma2 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma2<<G4endl;
G4cout<<"dsigma "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<dSigmadEprod<<G4endl;
}
//correction of differential cross section at high energy to correct for the suppression of particle at secondary at high
//energy used in the Bethe Bloch Model. This correction consist to multiply by g the probability function used
//to test the rejection of a secondary
//-------------------------
//Source code taken from G4BetheBlochModel::SampleSecondaries
G4double deltaKinEnergy = kinEnergyProd;
//Part of the taken code
//----------------------
// projectile formfactor - suppresion of high energy
// delta-electron production at high energy
G4double x = formfact*deltaKinEnergy;
if(x > 1.e-6) {
G4double totEnergy = kinEnergyProj + mass;
G4double etot2 = totEnergy*totEnergy;
G4double beta2 = kinEnergyProj*(kinEnergyProj + 2.0*mass)/etot2;
G4double f;
G4double f1 = 0.0;
f = 1.0 - beta2*deltaKinEnergy/Tmax;
if( 0.5 == spin ) {
f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
f += f1;
}
G4double x1 = 1.0 + x;
G4double g = 1.0/(x1*x1);
if( 0.5 == spin ) {
G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
g *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
}
if(g > 1.0) {
G4cout << "### G4BetheBlochModel in Adjoint Sim WARNING: g= " << g
<< G4endl;
g=1.;
}
//G4cout<<"g"<<g<<G4endl;
dSigmadEprod*=g;
}
}
return dSigmadEprod;
}
//////////////////////////////////////////////////////////////////////////////////////////////
//
void G4AdjointhIonisationModel::DefineProjectileProperty()
{
//Slightly modified code taken from G4BetheBlochModel::SetParticle
//------------------------------------------------
G4String pname = theDirectPrimaryPartDef->GetParticleName();
if (theDirectPrimaryPartDef->GetParticleType() == "nucleus" &&
pname != "deuteron" && pname != "triton") {
isIon = true;
}
mass = theDirectPrimaryPartDef->GetPDGMass();
spin = theDirectPrimaryPartDef->GetPDGSpin();
G4double q = theDirectPrimaryPartDef->GetPDGCharge()/eplus;
chargeSquare = q*q;
ratio = electron_mass_c2/mass;
ratio2 = ratio*ratio;
one_plus_ratio_2=(1+ratio)*(1+ratio);
one_minus_ratio_2=(1-ratio)*(1-ratio);
G4double magmom = theDirectPrimaryPartDef->GetPDGMagneticMoment()
*mass/(0.5*eplus*hbar_Planck*c_squared);
magMoment2 = magmom*magmom - 1.0;
formfact = 0.0;
if(theDirectPrimaryPartDef->GetLeptonNumber() == 0) {
G4double x = 0.8426*GeV;
if(spin == 0.0 && mass < GeV) {x = 0.736*GeV;}
else if(mass > GeV) {
x /= G4NistManager::Instance()->GetZ13(mass/proton_mass_c2);
// tlimit = 51.2*GeV*A13[iz]*A13[iz];
}
formfact = 2.0*electron_mass_c2/(x*x);
tlimit = 2.0/formfact;
}
}
////////////////////////////////////////////////////////////////////////////////
//
G4double G4AdjointhIonisationModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
G4double primEnergy,
G4bool IsScatProjToProjCase)
{
if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase);
DefineCurrentMaterial(aCouple);
G4double Cross=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2*mass;
if (!IsScatProjToProjCase ){
G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy);
G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy);
if (Emax_proj>Emin_proj && primEnergy > currentTcutForDirectSecond) {
Cross*=(1./Emin_proj -1./Emax_proj)/primEnergy;
}
else Cross=0.;
}
else {
G4double Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy);
G4double Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,currentTcutForDirectSecond);
G4double diff1=Emin_proj-primEnergy;
G4double diff2=Emax_proj-primEnergy;
G4double t1=(1./diff1+1./Emin_proj-1./diff2-1./Emax_proj)/primEnergy;
G4double t2=2.*std::log(diff2*Emin_proj/Emax_proj/diff1)/primEnergy/primEnergy;
Cross*=(t1+t2);
}
lastCS =Cross;
return Cross;
}
//////////////////////////////////////////////////////////////////////////////
//
G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
{
G4double Tmax=PrimAdjEnergy*one_plus_ratio_2/(one_minus_ratio_2-2.*ratio*PrimAdjEnergy/mass);
return Tmax;
}
//////////////////////////////////////////////////////////////////////////////
//
G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy,G4double Tcut)
{ return PrimAdjEnergy+Tcut;
}
//////////////////////////////////////////////////////////////////////////////
//
G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMaxForProdToProjCase(G4double )
{ return HighEnergyLimit;
}
//////////////////////////////////////////////////////////////////////////////
//
G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
{ G4double Tmin= (2*PrimAdjEnergy-4*mass + std::sqrt(4.*PrimAdjEnergy*PrimAdjEnergy +16.*mass*mass + 8.*PrimAdjEnergy*mass*(1/ratio +ratio)))/4.;
return Tmin;
}
source/processes/electromagnetic/adjoint/src/G4AdjointhMultipleScattering.cc
+149
View File
@@ -0,0 +1,149 @@
//
// ********************************************************************
// * License and Disclaimer *
// * *
// * The Geant4 software is copyright of the Copyright Holders of *
// * the Geant4 Collaboration. It is provided under the terms and *
// * conditions of the Geant4 Software License, included in the file *
// * LICENSE and available at http://cern.ch/geant4/license . These *
// * include a list of copyright holders. *
// * *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work make any representation or warranty, express or implied, *
// * regarding this software system or assume any liability for its *
// * use. Please see the license in the file LICENSE and URL above *
// * for the full disclaimer and the limitation of liability. *
// * *
// * This code implementation is the result of the scientific and *
// * technical work of the GEANT4 collaboration. *
// * By using, copying, modifying or distributing the software (or *
// * any work based on the software) you agree to acknowledge its *
// * use in resulting scientific publications, and indicate your *
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
// $Id: G4AdjointhMultipleScattering.cc,v 1.2 2009/11/20 10:31:20 ldesorgh Exp $
// GEANT4 tag $Name: geant4-09-03 $
//
//
// GEANT4 Class file
//
// File name: G4AdjointhMultipleScattering
//
// Author: Desorgher Laurent
//
// Creation date: 03.06.2009 cloned from G4hMultipleScattering by U.Laszlo with slight modification for adjoint_ion.
//
// -----------------------------------------------------------------------------
//
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#include "G4AdjointhMultipleScattering.hh"
#include "G4UrbanMscModel.hh"
#include "G4UrbanMscModel90.hh"
#include "G4MscStepLimitType.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
using namespace std;
G4AdjointhMultipleScattering::G4AdjointhMultipleScattering(const G4String& processName)
: G4VMultipleScattering(processName)
{
isInitialized = false;
isIon = false;
SetStepLimitType(fMinimal);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4AdjointhMultipleScattering::~G4AdjointhMultipleScattering()
{}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4bool G4AdjointhMultipleScattering::IsApplicable (const G4ParticleDefinition& p)
{
return (p.GetPDGCharge() != 0.0 && !p.IsShortLived());
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void G4AdjointhMultipleScattering::InitialiseProcess(const G4ParticleDefinition* p)
{
// Modification of parameters between runs
if(isInitialized) {
if (p->GetParticleType() != "adjoint_nucleus" && p->GetPDGMass() < GeV) {
mscUrban->SetStepLimitType(StepLimitType());
mscUrban->SetLateralDisplasmentFlag(LateralDisplasmentFlag());
mscUrban->SetSkin(Skin());
mscUrban->SetRangeFactor(RangeFactor());
mscUrban->SetGeomFactor(GeomFactor());
}
return;
}
// defaults for ions, which cannot be overwritten
if (p->GetParticleType() == "adjoint_nucleus" || p->GetPDGMass() > GeV) {
SetStepLimitType(fMinimal);
SetLateralDisplasmentFlag(false);
SetBuildLambdaTable(false);
if(p->GetParticleType() == "adjoint_nucleus") isIon = true;
}
// initialisation of parameters
G4String part_name = p->GetParticleName();
mscUrban = new G4UrbanMscModel90();
mscUrban->SetStepLimitType(StepLimitType());
mscUrban->SetLateralDisplasmentFlag(LateralDisplasmentFlag());
mscUrban->SetSkin(Skin());
mscUrban->SetRangeFactor(RangeFactor());
mscUrban->SetGeomFactor(GeomFactor());
AddEmModel(1,mscUrban);
isInitialized = true;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void G4AdjointhMultipleScattering::PrintInfo()
{
G4cout << " RangeFactor= " << RangeFactor()
<< ", step limit type: " << StepLimitType()
<< ", lateralDisplacement: " << LateralDisplasmentFlag()
<< ", skin= " << Skin()
// << ", geomFactor= " << GeomFactor()
<< G4endl;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
/*G4double G4AdjointhMultipleScattering::AlongStepGetPhysicalInteractionLength(
const G4Track& track,
double,
G4double currentMinimalStep,
G4double& currentSafety,
G4GPILSelection* selection)
{
// get Step limit proposed by the process
valueGPILSelectionMSC = NotCandidateForSelection;
G4double escaled = track.GetKineticEnergy();
if(isIon) escaled *= track.GetDynamicParticle()->GetMass()/proton_mass_c2;
G4double steplength = GetMscContinuousStepLimit(track,
escaled,
currentMinimalStep,
currentSafety);
// G4cout << "StepLimit= " << steplength << G4endl;
// set return value for G4GPILSelection
*selection = valueGPILSelectionMSC;
return steplength;
}
*/
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
source/processes/electromagnetic/adjoint/src/G4ContinuousGainOfEnergy.cc
+170 -28
View File
@@ -23,6 +23,9 @@
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
// $Id: G4ContinuousGainOfEnergy.cc,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
// GEANT4 tag $Name: geant4-09-03 $
//
#include "G4ContinuousGainOfEnergy.hh"
#include "G4Step.hh"
#include "G4ParticleDefinition.hh"
@@ -30,6 +33,8 @@
#include "G4VEmFluctuationModel.hh"
#include "G4VParticleChange.hh"
#include "G4UnitsTable.hh"
#include "G4AdjointCSManager.hh"
#include "G4LossTableManager.hh"
///////////////////////////////////////////////////////
@@ -38,11 +43,22 @@ G4ContinuousGainOfEnergy::G4ContinuousGainOfEnergy(const G4String& name,
G4ProcessType type): G4VContinuousProcess(name, type)
{
linLossLimit=0.05;
lossFluctuationArePossible =true;
lossFluctuationFlag=true;
is_integral = false;
//Will be properly set in SetDirectParticle()
IsIon=false;
massRatio =1.;
chargeSqRatio=1.;
preStepChargeSqRatio=1.;
}
///////////////////////////////////////////////////////
@@ -69,8 +85,20 @@ void G4ContinuousGainOfEnergy::BuildPhysicsTable(const G4ParticleDefinition&)
;
}
///////////////////////////////////////////////////////
//
void G4ContinuousGainOfEnergy::SetDirectParticle(G4ParticleDefinition* p)
{theDirectPartDef=p;
if (theDirectPartDef->GetParticleType()== "nucleus") {
IsIon=true;
massRatio = proton_mass_c2/theDirectPartDef->GetPDGMass();
G4double q=theDirectPartDef->GetPDGCharge();
chargeSqRatio=q*q;
}
}
///////////////////////////////////////////////////////
//
@@ -79,81 +107,131 @@ G4VParticleChange* G4ContinuousGainOfEnergy::AlongStepDoIt(const G4Track& track,
const G4Step& step)
{
//Caution in this method the step length should be the true step length
// A problem is that this is compute by the multiple scattering that does not know the energy at the end of the adjoint step. This energy is used during the
//Forward sim. Nothing we can really do against that at this time. This is inherent to the MS method
//
aParticleChange.Initialize(track);
// Get the actual (true) Step length
//----------------------------------
G4double length = step.GetStepLength();
G4double degain = 0.0;
// Compute this for weight change after continuous energy loss
//-------------------------------------------------------------
G4double DEDX_before =
theDirectEnergyLossProcess
->GetDEDX(preStepKinEnergy, currentCouple);
G4double DEDX_before = theDirectEnergyLossProcess->GetDEDX(preStepKinEnergy, currentCouple);
// For the fluctuation we generate a new dynamic particle with energy =preEnergy+egain
// and then compute the fluctuation given in the direct case.
//-----------------------------------------------------------------------
G4DynamicParticle* dynParticle = new G4DynamicParticle();
*dynParticle = *(track.GetDynamicParticle());
G4double Tkin = dynParticle->GetKineticEnergy();
dynParticle->SetDefinition(theDirectPartDef);
G4double Tkin = dynParticle->GetKineticEnergy();
G4double Tkin1=Tkin*0.001;
size_t n=1;
if (is_integral ) n=10;
n=1;
G4double dlength= length/n;
for (size_t i=0;i<n;i++) {
G4double factor_dE=1.;
if (Tkin != preStepKinEnergy && IsIon) {
chargeSqRatio = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,Tkin);
theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatio);
}
G4double r = theDirectEnergyLossProcess->GetRange(Tkin, currentCouple);
if( dlength <= linLossLimit * r ) {
degain = DEDX_before*dlength;
G4double degain1 = dlength*theDirectEnergyLossProcess->GetDEDX(Tkin1, currentCouple);
factor_dE=1.+(degain1-degain)/(Tkin1-Tkin);
}
else {
G4double x = r + length;
degain = theDirectEnergyLossProcess->GetKineticEnergy(x,currentCouple) - theDirectEnergyLossProcess->GetKineticEnergy(r,currentCouple);
G4double x = r + dlength;
//degain = theDirectEnergyLossProcess->GetKineticEnergy(x,currentCouple) - theDirectEnergyLossProcess->GetKineticEnergy(r,currentCouple);
G4double E = theDirectEnergyLossProcess->GetKineticEnergy(x,currentCouple);
if (IsIon){
chargeSqRatio = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,E);
theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatio);
G4double x1= theDirectEnergyLossProcess->GetRange(E, currentCouple);
while (std::abs(x-x1)>0.01*x) {
E = theDirectEnergyLossProcess->GetKineticEnergy(x,currentCouple);
chargeSqRatio = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,E);
theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatio);
x1= theDirectEnergyLossProcess->GetRange(E, currentCouple);
}
}
G4double r1 = theDirectEnergyLossProcess->GetRange(Tkin1, currentCouple);
G4double x1 = r1 + dlength;
G4double E1 = theDirectEnergyLossProcess->GetKineticEnergy(x1,currentCouple);
factor_dE=(E1-E)/(Tkin1-Tkin);
degain=E-Tkin;
}
G4VEmModel* currentModel = theDirectEnergyLossProcess->SelectModelForMaterial(Tkin+degain,currentMaterialIndex);
//G4cout<<degain<<G4endl;
G4double tmax = currentModel->MaxSecondaryKinEnergy(dynParticle);
tmax = std::min(tmax,currentTcut);
dynParticle->SetKineticEnergy(Tkin+degain);
// Corrections, which cannot be tabulated for ions
//----------------------------------------
G4double esecdep=0;//not used in most models
currentModel->CorrectionsAlongStep(currentCouple, dynParticle, degain,esecdep, dlength);
// Sample fluctuations
//-------------------
G4double deltaE =0.;
if (lossFluctuationFlag ) {
deltaE = currentModel->GetModelOfFluctuations()->
SampleFluctuations(currentMaterial,dynParticle,tmax,length,degain)-degain;
SampleFluctuations(currentMaterial,dynParticle,tmax,dlength,degain)-degain;
}
Tkin+=degain+deltaE;
G4double egain=degain+deltaE;
if (egain <=0) egain=degain;
Tkin+=egain;
dynParticle->SetKineticEnergy(Tkin);
}
}
// Corrections, which cannot be tabulated
// probably this should be also changed
// at this time it does nothing so we can leave it
//CorrectionsAlongStep(currentCouple, dynParticle, egain, length);
delete dynParticle;
if (IsIon){
chargeSqRatio = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,Tkin);
theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatio);
}
G4double DEDX_after = theDirectEnergyLossProcess->GetDEDX(Tkin, currentCouple);
G4double weight_correction=DEDX_after/DEDX_before; //probably not needed
weight_correction=1.;
G4double weight_correction=DEDX_after/DEDX_before;
aParticleChange.ProposeEnergy(Tkin);
//we still need to register in the particleChange the modification of the weight of the particle
G4double new_weight=weight_correction*track.GetWeight();
aParticleChange.SetParentWeightByProcess(true);
aParticleChange.SetParentWeightByProcess(false);
aParticleChange.ProposeParentWeight(new_weight);
@@ -167,3 +245,67 @@ void G4ContinuousGainOfEnergy::SetLossFluctuations(G4bool val)
if(val && !lossFluctuationArePossible) return;
lossFluctuationFlag = val;
}
///////////////////////////////////////////////////////
//
G4double G4ContinuousGainOfEnergy::GetContinuousStepLimit(const G4Track& track,
G4double , G4double , G4double& )
{
G4double x = DBL_MAX;
x=.1*mm;
DefineMaterial(track.GetMaterialCutsCouple());
preStepKinEnergy = track.GetKineticEnergy();
preStepScaledKinEnergy = track.GetKineticEnergy()*massRatio;
currentModel = theDirectEnergyLossProcess->SelectModelForMaterial(preStepScaledKinEnergy,currentCoupleIndex);
G4double emax_model=currentModel->HighEnergyLimit();
if (IsIon) {
chargeSqRatio = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,preStepKinEnergy);
preStepChargeSqRatio = chargeSqRatio;
theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,preStepChargeSqRatio);
}
G4double maxE =1.1*preStepKinEnergy;
/*if (preStepKinEnergy< 0.05*MeV) maxE =2.*preStepKinEnergy;
else if (preStepKinEnergy< 0.1*MeV) maxE =1.5*preStepKinEnergy;
else if (preStepKinEnergy< 0.5*MeV) maxE =1.25*preStepKinEnergy;*/
if (preStepKinEnergy < currentTcut) maxE = std::min(currentTcut,maxE);
maxE=std::min(emax_model*1.001,maxE);
G4double r = theDirectEnergyLossProcess->GetRange(preStepKinEnergy, currentCouple);
if (IsIon) {
G4double chargeSqRatioAtEmax = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,maxE);
theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatioAtEmax);
}
G4double r1 = theDirectEnergyLossProcess->GetRange(maxE, currentCouple);
if (IsIon) theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,preStepChargeSqRatio);
x=r1-r;
x=std::max(r1-r,0.001*mm);
return x;
}
#include "G4EmCorrections.hh"
///////////////////////////////////////////////////////
//
void G4ContinuousGainOfEnergy::SetDynamicMassCharge(const G4Track& ,G4double energy)
{
G4double ChargeSqRatio= G4LossTableManager::Instance()->EmCorrections()->EffectiveChargeSquareRatio(theDirectPartDef,currentMaterial,energy);
if (theDirectEnergyLossProcess) theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,ChargeSqRatio);
}
source/processes/electromagnetic/adjoint/src/G4InversePEEffect.cc
+8 -2
View File
@@ -23,6 +23,9 @@
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
// $Id: G4InversePEEffect.cc,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
// GEANT4 tag $Name: geant4-09-03 $
//
#include "G4InversePEEffect.hh"
#include "G4VEmAdjointModel.hh"
#include "G4AdjointPhotoElectricModel.hh"
@@ -30,9 +33,12 @@
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
G4InversePEEffect::G4InversePEEffect(G4String process_name,G4AdjointPhotoElectricModel* aModel):
G4VAdjointInverseScattering(process_name,false)
G4VAdjointReverseReaction(process_name,false)
{theAdjointEMModel = aModel;
theAdjointEMModel->SetSecondPartOfSameType(false);
theAdjointEMModel->SetSecondPartOfSameType(false);
SetIntegralMode(false);
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
source/processes/electromagnetic/adjoint/src/G4IonInverseIonisation.cc
+51
View File
@@ -0,0 +1,51 @@
//
// ********************************************************************
// * License and Disclaimer *
// * *
// * The Geant4 software is copyright of the Copyright Holders of *
// * the Geant4 Collaboration. It is provided under the terms and *
// * conditions of the Geant4 Software License, included in the file *
// * LICENSE and available at http://cern.ch/geant4/license . These *
// * include a list of copyright holders. *
// * *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work make any representation or warranty, express or implied, *
// * regarding this software system or assume any liability for its *
// * use. Please see the license in the file LICENSE and URL above *
// * for the full disclaimer and the limitation of liability. *
// * *
// * This code implementation is the result of the scientific and *
// * technical work of the GEANT4 collaboration. *
// * By using, copying, modifying or distributing the software (or *
// * any work based on the software) you agree to acknowledge its *
// * use in resulting scientific publications, and indicate your *
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
// $Id: G4IonInverseIonisation.cc,v 1.2 2009/11/20 10:31:20 ldesorgh Exp $
// GEANT4 tag $Name: geant4-09-03 $
//
///////////////////////////////////////////////////////
// File name: G4IonInverseIonisation
//
// Author: Laurent Desorgher
//
// Creation date: 25.08.2009
//
///////////////////////////////////////////////////////
#include "G4IonInverseIonisation.hh"
#include "G4VEmAdjointModel.hh"
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
G4IonInverseIonisation::G4IonInverseIonisation(G4bool whichScatCase,G4String process_name,G4AdjointIonIonisationModel* aEmAdjointModel):
G4VAdjointReverseReaction(process_name,whichScatCase)
{theAdjointEMModel = aEmAdjointModel;
theAdjointEMModel->SetSecondPartOfSameType(false);
SetIntegralMode(true);
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
G4IonInverseIonisation::~G4IonInverseIonisation(){
}
source/processes/electromagnetic/adjoint/src/G4VAdjointInverseScattering.cc → source/processes/electromagnetic/adjoint/src/G4VAdjointReverseReaction.cc
+48 -13
View File
@@ -23,7 +23,10 @@
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
#include "G4VAdjointInverseScattering.hh"
// $Id: G4VAdjointReverseReaction.cc,v 1.2 2009/11/20 10:31:20 ldesorgh Exp $
// GEANT4 tag $Name: geant4-09-03 $
//
#include "G4VAdjointReverseReaction.hh"
#include "G4AdjointCSManager.hh"
#include "G4AdjointCSMatrix.hh"
#include "G4AdjointInterpolator.hh"
@@ -38,30 +41,28 @@
#include "G4AdjointElectron.hh"
G4VAdjointInverseScattering::
G4VAdjointInverseScattering(G4String process_name, G4bool whichScatCase):
G4VAdjointReverseReaction::
G4VAdjointReverseReaction(G4String process_name, G4bool whichScatCase):
G4VDiscreteProcess(process_name)
{theAdjointCSManager = G4AdjointCSManager::GetAdjointCSManager();
IsScatProjToProjCase=whichScatCase;
/*theAdjointEMModel=aModel;
IsScatProjToProjCase=whichScatCase;*/
fParticleChange=new G4ParticleChange();
}
//////////////////////////////////////////////////////////////////////////////
//
G4VAdjointInverseScattering::
~G4VAdjointInverseScattering()
G4VAdjointReverseReaction::
~G4VAdjointReverseReaction()
{;
}
//////////////////////////////////////////////////////////////////////////////
//
void G4VAdjointInverseScattering::PreparePhysicsTable(const G4ParticleDefinition&)
void G4VAdjointReverseReaction::PreparePhysicsTable(const G4ParticleDefinition&)
{;
}
//////////////////////////////////////////////////////////////////////////////
//
void G4VAdjointInverseScattering::BuildPhysicsTable(const G4ParticleDefinition&)
void G4VAdjointReverseReaction::BuildPhysicsTable(const G4ParticleDefinition&)
{
theAdjointCSManager->BuildCrossSectionMatrices(); //do not worry it will be done just once
@@ -70,10 +71,25 @@ void G4VAdjointInverseScattering::BuildPhysicsTable(const G4ParticleDefinition&)
}
//////////////////////////////////////////////////////////////////////////////
//
G4VParticleChange* G4VAdjointInverseScattering::PostStepDoIt(const G4Track& track, const G4Step& )
G4VParticleChange* G4VAdjointReverseReaction::PostStepDoIt(const G4Track& track, const G4Step& )
{
fParticleChange->Initialize(track);
/* if (IsFwdCSUsed && IsIntegralModeUsed){ //INtegral mode still unstable
G4double Tkin = step.GetPostStepPoint()->GetKineticEnergy();
G4double fwdCS = theAdjointCSManager->GetTotalForwardCS(track.GetDefinition(), Tkin, track.GetMaterialCutsCouple());
//G4cout<<"lastCS "<<lastCS<<G4endl;
if (fwdCS<lastCS*G4UniformRand()) { // the reaction does not take place, same integral method as the one used for forward ionisation in G4
ClearNumberOfInteractionLengthLeft();
return fParticleChange;
}
}
*/
theAdjointEMModel->SampleSecondaries(track,
IsScatProjToProjCase,
fParticleChange);
@@ -86,7 +102,7 @@ G4VParticleChange* G4VAdjointInverseScattering::PostStepDoIt(const G4Track& trac
}
//////////////////////////////////////////////////////////////////////////////
//
G4double G4VAdjointInverseScattering::GetMeanFreePath(const G4Track& track,
G4double G4VAdjointReverseReaction::GetMeanFreePath(const G4Track& track,
G4double ,
G4ForceCondition* condition)
{ *condition = NotForced;
@@ -95,9 +111,28 @@ G4double G4VAdjointInverseScattering::GetMeanFreePath(const G4Track& track,
G4double Sigma =
theAdjointEMModel->AdjointCrossSection(track.GetMaterialCutsCouple(),preStepKinEnergy,IsScatProjToProjCase);
G4double fwd_TotCS;
Sigma *= theAdjointCSManager->GetCrossSectionCorrection(track.GetDefinition(),preStepKinEnergy,track.GetMaterialCutsCouple(),IsFwdCSUsed, fwd_TotCS);
//G4cout<<fwd_TotCS<<G4endl;
/*if (IsFwdCSUsed && IsIntegralModeUsed){ //take the maximum cross section only for charged particle
G4double e_sigma_max, sigma_max;
theAdjointCSManager->GetMaxFwdTotalCS(track.GetDefinition(),
track.GetMaterialCutsCouple(), e_sigma_max, sigma_max);
if (e_sigma_max > preStepKinEnergy){
Sigma*=sigma_max/fwd_TotCS;
}
}
*/
G4double mean_free_path = 1.e60 *mm;
if (Sigma>0) mean_free_path = 1./Sigma;
lastCS=Sigma;
/*G4cout<<"Sigma "<<Sigma<<G4endl;
G4cout<<"mean_free_path [mm] "<<mean_free_path/mm<<G4endl;
*/
G4double mean_free_path = 1./Sigma;
//G4cout<<"mean_free_path [mm] "<<mean_free_path/mm<<std::endl;
return mean_free_path;
}
source/processes/electromagnetic/adjoint/src/G4VEmAdjointModel.cc
+210 -377
View File
@@ -23,213 +23,60 @@
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
// $Id: G4VEmAdjointModel.cc,v 1.5 2009/12/16 17:50:09 gunter Exp $
// GEANT4 tag $Name: geant4-09-03 $
//
#include "G4VEmAdjointModel.hh"
#include "G4AdjointCSManager.hh"
#include "G4Integrator.hh"
#include "G4TrackStatus.hh"
#include "G4ParticleChange.hh"
#include "G4AdjointElectron.hh"
#include "G4AdjointInterpolator.hh"
#include "G4PhysicsTable.hh"
////////////////////////////////////////////////////////////////////////////////
//
G4VEmAdjointModel::G4VEmAdjointModel(const G4String& nam):
name(nam)
// lowLimit(0.1*keV), highLimit(100.0*TeV), fluc(0), name(nam), pParticleChange(0)
{ G4AdjointCSManager::GetAdjointCSManager()->RegisterEmAdjointModel(this);
CorrectWeightMode =true;
UseMatrix =true;
UseMatrixPerElement = true;
ApplyCutInRange = true;
ApplyBiasing = true;
UseOnlyOneMatrixForAllElements = true;
IsIonisation =true;
CS_biasing_factor =1.;
//ApplyBiasing = false;
{
G4AdjointCSManager::GetAdjointCSManager()->RegisterEmAdjointModel(this);
second_part_of_same_type =false;
theDirectEMModel=0;
}
////////////////////////////////////////////////////////////////////////////////
//
G4VEmAdjointModel::~G4VEmAdjointModel()
{;}
////////////////////////////////////////////////////////////////////////////////
//
void G4VEmAdjointModel::SampleSecondaries(const G4Track& aTrack,
G4bool IsScatProjToProjCase,
G4ParticleChange* fParticleChange)
{
const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
//DefineCurrentMaterial(aTrack->GetMaterialCutsCouple());
size_t ind=0;
if (!UseMatrixPerElement) ind = currentMaterialIndex;
//G4cout<<theAdjointPrimary<<std::endl;
else if (!UseOnlyOneMatrixForAllElements) { //Select Material
std::vector<double>* CS_Vs_Element = &CS_Vs_ElementForScatProjToProjCase;
if ( !IsScatProjToProjCase) CS_Vs_Element = &CS_Vs_ElementForProdToProjCase;
G4double rand_var= G4UniformRand();
G4double SumCS=0.;
for (size_t i=0;i<CS_Vs_Element->size();i++){
SumCS+=(*CS_Vs_Element)[i];
if (rand_var<=SumCS/lastCS){
ind=i;
break;
}
}
ind = currentMaterial->GetElement(ind)->GetIndex();
}
//Elastic inverse scattering //not correct in all the cases
//---------------------------------------------------------
G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy();
G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
//G4cout<<adjointPrimKinEnergy<<std::endl;
if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
return;
}
//Sample secondary energy
//-----------------------
G4double projectileKinEnergy;
// if (!IsIonisation ) {
projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(ind,
adjointPrimKinEnergy,
IsScatProjToProjCase);
//}
/*else {
projectileKinEnergy = SampleAdjSecEnergyFromDiffCrossSectionPerAtom(adjointPrimKinEnergy,IsScatProjToProjCase);
//G4cout<<projectileKinEnergy<<std::endl;
}*/
//Weight correction
//-----------------------
CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(), adjointPrimKinEnergy,projectileKinEnergy);
//Kinematic
//---------
G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
//Companion
//-----------
G4double companionM0;
companionM0=(adjointPrimTotalEnergy-adjointPrimKinEnergy);
if (IsScatProjToProjCase) {
companionM0=theAdjEquivOfDirectSecondPartDef->GetPDGMass();
}
G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;
//Projectile momentum
//--------------------
G4double P_parallel = (adjointPrimP*adjointPrimP + projectileP2 - companionP2)/(2.*adjointPrimP);
G4double P_perp = std::sqrt( projectileP2 - P_parallel*P_parallel);
G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
G4double phi =G4UniformRand()*2.*3.1415926;
G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
projectileMomentum.rotateUz(dir_parallel);
if (!IsScatProjToProjCase && CorrectWeightMode){ //kill the primary and add a secondary
fParticleChange->ProposeTrackStatus(fStopAndKill);
fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
//G4cout<<"projectileMomentum "<<projectileMomentum<<std::endl;
}
else {
fParticleChange->ProposeEnergy(projectileKinEnergy);
fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
}
}
////////////////////////////////////////////////////////////////////////////////
//
void G4VEmAdjointModel::CorrectPostStepWeight(G4ParticleChange* fParticleChange, G4double old_weight, G4double , G4double )
{
G4double new_weight=old_weight;
if (CorrectWeightMode) {
G4double w_corr =1./CS_biasing_factor;
//G4cout<<w_corr<<std::endl;
/*G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection(theAdjEquivOfDirectPrimPartDef,
theAdjEquivOfDirectSecondPartDef,
adjointPrimKinEnergy,projectileKinEnergy,
aTrack.GetMaterialCutsCouple());
w_corr = projectileKinEnergy;
G4double Emin,Emax;
if (IsScatProjToProjCase) {
Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy, currentTcutForDirectSecond);
}
else {
Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
Emin = GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);
}
w_corr *=std::log(Emax/Emin)/(Emax-Emin); */
new_weight*=w_corr;
}
G4cout<< "new weight"<<new_weight<<std::endl;
fParticleChange->SetParentWeightByProcess(false);
fParticleChange->SetSecondaryWeightByProcess(false);
fParticleChange->ProposeParentWeight(new_weight);
}
////////////////////////////////////////////////////////////////////////////////
//
G4double G4VEmAdjointModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
G4double primEnergy,
G4bool IsScatProjToProjCase)
{
DefineCurrentMaterial(aCouple);
//G4double fwdCS = G4AdjointCSManager::GetAdjointCSManager()->GetTotalForwardCS(G4AdjointElectron::AdjointElectron(),primEnergy,aCouple);
//G4double adjCS = G4AdjointCSManager::GetAdjointCSManager()->GetTotalAdjointCS(G4AdjointElectron::AdjointElectron(), primEnergy,aCouple);
if (IsScatProjToProjCase){
lastCS = G4AdjointCSManager::GetAdjointCSManager()->ComputeAdjointCS(currentMaterial,
this,
primEnergy,
currentTcutForDirectSecond,
true,
CS_Vs_ElementForScatProjToProjCase);
/*G4double fwdCS = G4AdjointCSManager::GetAdjointCSManager()->GetTotalForwardCS(theAdjEquivOfDirectPrimPartDef,primEnergy,aCouple);
G4double adjCS = G4AdjointCSManager::GetAdjointCSManager()->GetTotalAdjointCS(theAdjEquivOfDirectPrimPartDef, primEnergy,aCouple);
*/
//if (adjCS >0 )lastCS *=fwdCS/adjCS;
}
else {
lastCS = G4AdjointCSManager::GetAdjointCSManager()->ComputeAdjointCS(currentMaterial,
this,
primEnergy,
currentTcutForDirectSecond,
false,
CS_Vs_ElementForProdToProjCase);
/*G4double fwdCS = G4AdjointCSManager::GetAdjointCSManager()->GetTotalForwardCS(theAdjEquivOfDirectSecondPartDef,primEnergy,aCouple);
G4double adjCS = G4AdjointCSManager::GetAdjointCSManager()->GetTotalAdjointCS(theAdjEquivOfDirectSecondPartDef, primEnergy,aCouple);
*/
//if (adjCS >0 )lastCS *=fwdCS/adjCS;
//lastCS=0.;
}
preStepEnergy=primEnergy;
std::vector<G4double>* CS_Vs_Element = &CS_Vs_ElementForProdToProjCase;
if (IsScatProjToProjCase) CS_Vs_Element = &CS_Vs_ElementForScatProjToProjCase;
lastCS = G4AdjointCSManager::GetAdjointCSManager()->ComputeAdjointCS(currentMaterial,
this,
primEnergy,
currentTcutForDirectSecond,
IsScatProjToProjCase,
*CS_Vs_Element);
if (IsScatProjToProjCase) lastAdjointCSForScatProjToProjCase = lastCS;
else lastAdjointCSForProdToProjCase =lastCS;
return lastCS;
}
////////////////////////////////////////////////////////////////////////////////
//
//The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
//General implementation correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
G4double G4VEmAdjointModel::DiffCrossSectionPerAtomPrimToSecond(
G4double kinEnergyProj,
G4double kinEnergyProd,
@@ -244,31 +91,15 @@ G4double G4VEmAdjointModel::DiffCrossSectionPerAtomPrimToSecond(
if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ //the produced particle should have a kinetic energy smaller than the projectile
G4double Tmax=kinEnergyProj;
if (second_part_of_same_type) Tmax = kinEnergyProj/2.;
return Z*DiffCrossSectionMoller(kinEnergyProj,kinEnergyProd);
//it could be thta Tmax here should be DBLMAX
//Tmax=DBLMAX;
G4double E1=kinEnergyProd;
G4double E1=kinEnergyProd;
G4double E2=kinEnergyProd*1.000001;
G4double dE=(E2-E1);
G4double sigma1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
G4double sigma2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
dSigmadEprod=(sigma1-sigma2)/dE;
if (dSigmadEprod>1.) {
G4cout<<"sigma1 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma1<<std::endl;
G4cout<<"sigma2 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma2<<std::endl;
G4cout<<"dsigma "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<dSigmadEprod<<std::endl;
}
}
return dSigmadEprod;
@@ -306,33 +137,13 @@ G4double G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToSecond(
if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){
G4double Tmax=kinEnergyProj;
if (second_part_of_same_type) Tmax = kinEnergyProj/2.;
//it could be thta Tmax here should be DBLMAX
//Tmax=DBLMAX;
G4double E1=kinEnergyProd;
G4double E2=kinEnergyProd*1.0001;
G4double E1=kinEnergyProd;
G4double E2=kinEnergyProd*1.0001;
G4double dE=(E2-E1);
G4double sigma1=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E1,E2);
//G4double sigma2=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E2,1.e50);
dSigmadEprod=sigma1/dE;
if (dSigmadEprod <0) { //could happen with bremstrahlung dur to suppression effect
G4cout<<"Halllllllllllllllllllllllllllllllllllllllllllllllo "<<kinEnergyProj<<'\t'<<E1<<'\t'<<dSigmadEprod<<std::endl;
E1=kinEnergyProd;
E2=E1*1.1;
dE=E2-E1;
sigma1=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E1,1.e50);
G4double sigma2=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E2,1.e50);
dSigmadEprod=(sigma1-sigma2)/dE;
G4cout<<dSigmadEprod<<std::endl;
}
G4double sigma2=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E2,1.e50);
dSigmadEprod=(sigma1-sigma2)/dE;
}
return dSigmadEprod;
@@ -355,49 +166,28 @@ G4double G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToScatPrim(
///////////////////////////////////////////////////////////////////////////////////////////////////////////
//
G4double G4VEmAdjointModel::DiffCrossSectionFunction1(G4double kinEnergyProj){
//return kinEnergyProj*kinEnergyProj;
//ApplyBiasing=false;
G4double bias_factor = CS_biasing_factor*kinEnergyProdForIntegration/kinEnergyProj;
if (!ApplyBiasing) bias_factor =CS_biasing_factor;
//G4cout<<bias_factor<<std::endl;
if (UseMatrixPerElement ) {
return DiffCrossSectionPerAtomPrimToSecond(kinEnergyProj,kinEnergyProdForIntegration,ZSelectedNucleus,ASelectedNucleus)*bias_factor;
}
else {
else {
return DiffCrossSectionPerVolumePrimToSecond(SelectedMaterial,kinEnergyProj,kinEnergyProdForIntegration)*bias_factor;
}
}
//////////////////////////////////////////////////////////////////////////////
//
G4double G4VEmAdjointModel::DiffCrossSectionMoller(G4double kinEnergyProj,G4double kinEnergyProd){
G4double electron_mass_c2=0.51099906*MeV;
G4double energy = kinEnergyProj + electron_mass_c2;
G4double x = kinEnergyProd/kinEnergyProj;
G4double gam = energy/electron_mass_c2;
G4double gamma2 = gam*gam;
G4double beta2 = 1.0 - 1.0/gamma2;
G4double g = (2.0*gam - 1.0)/gamma2;
G4double y = 1.0 - x;
G4double fac=twopi_mc2_rcl2/electron_mass_c2;
G4double dCS = fac*( 1.-g + ((1.0 - g*x)/(x*x)) + ((1.0 - g*y)/(y*y)))/(beta2*(gam-1));
return dCS/kinEnergyProj;
}
////////////////////////////////////////////////////////////////////////////////
//
G4double G4VEmAdjointModel::DiffCrossSectionFunction2(G4double kinEnergyProj){
//return kinEnergyProj*kinEnergyProj;
G4double bias_factor = CS_biasing_factor*kinEnergyScatProjForIntegration/kinEnergyProj;
//ApplyBiasing=false;
if (!ApplyBiasing) bias_factor = CS_biasing_factor;
//G4cout<<bias_factor<<std::endl;
G4double bias_factor = CS_biasing_factor*kinEnergyScatProjForIntegration/kinEnergyProj;
if (UseMatrixPerElement ) {
return DiffCrossSectionPerAtomPrimToScatPrim(kinEnergyProj,kinEnergyScatProjForIntegration,ZSelectedNucleus,ASelectedNucleus)*bias_factor;
}
else {
else {
return DiffCrossSectionPerVolumePrimToScatPrim(SelectedMaterial,kinEnergyProj,kinEnergyScatProjForIntegration)*bias_factor;
}
@@ -405,14 +195,22 @@ G4double G4VEmAdjointModel::DiffCrossSectionFunction2(G4double kinEnergyProj){
}
////////////////////////////////////////////////////////////////////////////////
//
G4double G4VEmAdjointModel::DiffCrossSectionPerVolumeFunctionForIntegrationOverEkinProj(G4double kinEnergyProd)
{
return DiffCrossSectionPerVolumePrimToSecond(SelectedMaterial,kinEnergyProjForIntegration,kinEnergyProd);
}
////////////////////////////////////////////////////////////////////////////////
//
std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerAtomForSecond(
G4double kinEnergyProd,
G4double Z,
G4double A ,
G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
{ G4Integrator<G4VEmAdjointModel, G4double(G4VEmAdjointModel::*)(G4double)> integral;
ASelectedNucleus= G4int(A);
ZSelectedNucleus=G4int(Z);
{
G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
ASelectedNucleus= int(A);
ZSelectedNucleus=int(Z);
kinEnergyProdForIntegration = kinEnergyProd;
//compute the vector of integrated cross sections
@@ -421,24 +219,24 @@ std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSect
G4double minEProj= GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
G4double maxEProj= GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
G4double E1=minEProj;
std::vector< G4double >* log_ESec_vector = new std::vector< G4double >();
std::vector< G4double >* log_Prob_vector = new std::vector< G4double >();
std::vector< double>* log_ESec_vector = new std::vector< double>();
std::vector< double>* log_Prob_vector = new std::vector< double>();
log_ESec_vector->clear();
log_Prob_vector->clear();
log_ESec_vector->push_back(std::log(E1));
log_Prob_vector->push_back(-50.);
G4double E2=std::pow(10.,G4double( G4int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
G4double E2=std::pow(10.,double( int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
G4double fE=std::pow(10.,1./nbin_pro_decade);
G4double int_cross_section=0.;
if (std::pow(fE,5.)>(maxEProj/minEProj)) fE = std::pow(maxEProj/minEProj,0.2);
while (E1 <maxEProj*0.9999999){
//G4cout<<E1<<'\t'<<E2<<std::endl;
//G4cout<<E1<<'\t'<<E2<<G4endl;
int_cross_section +=integral.Simpson(this, &G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 10);
//G4cout<<"int_cross_section 1 "<<'\t'<<int_cross_section<<std::endl;
int_cross_section +=integral.Simpson(this,
&G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 5);
log_ESec_vector->push_back(std::log(std::min(E2,maxEProj)));
log_Prob_vector->push_back(std::log(int_cross_section));
E1=E2;
@@ -462,9 +260,9 @@ std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSect
G4double Z,
G4double A ,
G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
{ G4Integrator<G4VEmAdjointModel, G4double(G4VEmAdjointModel::*)(G4double)> integral;
ASelectedNucleus=G4int(A);
ZSelectedNucleus=G4int(Z);
{ G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
ASelectedNucleus=int(A);
ZSelectedNucleus=int(Z);
kinEnergyScatProjForIntegration = kinEnergyScatProj;
//compute the vector of integrated cross sections
@@ -478,29 +276,29 @@ std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSect
G4double dE2=dEmin;
std::vector< G4double >* log_ESec_vector = new std::vector< G4double >();
std::vector< G4double >* log_Prob_vector = new std::vector< G4double >();
std::vector< double>* log_ESec_vector = new std::vector< double>();
std::vector< double>* log_Prob_vector = new std::vector< double>();
log_ESec_vector->push_back(std::log(dEmin));
log_Prob_vector->push_back(-50.);
G4int nbins=std::max( G4int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
G4int nbins=std::max( int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
G4double fE=std::pow(dEmax/dEmin,1./nbins);
G4double int_cross_section=0.;
while (dE1 <dEmax*0.9999999999999){
dE2=dE1*fE;
int_cross_section +=integral.Simpson(this,
&G4VEmAdjointModel::DiffCrossSectionFunction2,minEProj+dE1,std::min(minEProj+dE2,maxEProj), 20);
//G4cout<<"int_cross_section "<<minEProj+dE1<<'\t'<<int_cross_section<<std::endl;
log_ESec_vector->push_back(std::log(std::min(dE2,maxEProj)));
&G4VEmAdjointModel::DiffCrossSectionFunction2,minEProj+dE1,std::min(minEProj+dE2,maxEProj), 5);
//G4cout<<"int_cross_section "<<minEProj+dE1<<'\t'<<int_cross_section<<G4endl;
log_ESec_vector->push_back(std::log(std::min(dE2,maxEProj-minEProj)));
log_Prob_vector->push_back(std::log(int_cross_section));
dE1=dE2;
}
/*G4cout<<"total int_cross_section"<<'\t'<<int_cross_section<<std::endl;
G4cout<<"energy "<<kinEnergyScatProj<<std::endl;*/
std::vector< std::vector<G4double> *> res_mat;
@@ -518,35 +316,32 @@ std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSect
G4Material* aMaterial,
G4double kinEnergyProd,
G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
{ G4Integrator<G4VEmAdjointModel, G4double(G4VEmAdjointModel::*)(G4double)> integral;
{ G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
SelectedMaterial= aMaterial;
kinEnergyProdForIntegration = kinEnergyProd;
//G4cout<<aMaterial->GetName()<<std::endl;
//G4cout<<kinEnergyProd/MeV<<std::endl;
//compute the vector of integrated cross sections
//compute the vector of integrated cross sections
//-------------------
G4double minEProj= GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
G4double maxEProj= GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
G4double E1=minEProj;
std::vector< G4double >* log_ESec_vector = new std::vector< G4double >();
std::vector< G4double >* log_Prob_vector = new std::vector< G4double >();
std::vector< double>* log_ESec_vector = new std::vector< double>();
std::vector< double>* log_Prob_vector = new std::vector< double>();
log_ESec_vector->clear();
log_Prob_vector->clear();
log_ESec_vector->push_back(std::log(E1));
log_Prob_vector->push_back(-50.);
G4double E2=std::pow(10.,G4double( G4int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
G4double E2=std::pow(10.,double( int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
G4double fE=std::pow(10.,1./nbin_pro_decade);
G4double int_cross_section=0.;
if (std::pow(fE,5.)>(maxEProj/minEProj)) fE = std::pow(maxEProj/minEProj,0.2);
while (E1 <maxEProj*0.9999999){
//G4cout<<E1<<'\t'<<E2<<std::endl;
int_cross_section +=integral.Simpson(this, &G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 10);
//G4cout<<"int_cross_section 1 "<<E1<<'\t'<<int_cross_section<<std::endl;
int_cross_section +=integral.Simpson(this,
&G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 5);
log_ESec_vector->push_back(std::log(std::min(E2,maxEProj)));
log_Prob_vector->push_back(std::log(int_cross_section));
E1=E2;
@@ -556,10 +351,12 @@ std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSect
std::vector< std::vector<G4double>* > res_mat;
res_mat.clear();
//if (int_cross_section >0.) {
if (int_cross_section >0.) {
res_mat.push_back(log_ESec_vector);
res_mat.push_back(log_Prob_vector);
//}
}
return res_mat;
}
@@ -570,12 +367,10 @@ std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSect
G4Material* aMaterial,
G4double kinEnergyScatProj,
G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
{ G4Integrator<G4VEmAdjointModel, G4double(G4VEmAdjointModel::*)(G4double)> integral;
{ G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
SelectedMaterial= aMaterial;
kinEnergyScatProjForIntegration = kinEnergyScatProj;
/*G4cout<<name<<std::endl;
G4cout<<aMaterial->GetName()<<std::endl;
G4cout<<kinEnergyScatProj/MeV<<std::endl;*/
//compute the vector of integrated cross sections
//-------------------
@@ -589,11 +384,11 @@ std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSect
G4double dE2=dEmin;
std::vector< G4double >* log_ESec_vector = new std::vector< G4double >();
std::vector< G4double >* log_Prob_vector = new std::vector< G4double >();
std::vector< double>* log_ESec_vector = new std::vector< double>();
std::vector< double>* log_Prob_vector = new std::vector< double>();
log_ESec_vector->push_back(std::log(dEmin));
log_Prob_vector->push_back(-50.);
G4int nbins=std::max( G4int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
G4int nbins=std::max( int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
G4double fE=std::pow(dEmax/dEmin,1./nbins);
G4double int_cross_section=0.;
@@ -601,9 +396,8 @@ std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSect
while (dE1 <dEmax*0.9999999999999){
dE2=dE1*fE;
int_cross_section +=integral.Simpson(this,
&G4VEmAdjointModel::DiffCrossSectionFunction2,minEProj+dE1,std::min(minEProj+dE2,maxEProj), 20);
//G4cout<<"int_cross_section "<<minEProj+dE1<<'\t'<<int_cross_section<<std::endl;
log_ESec_vector->push_back(std::log(std::min(dE2,maxEProj)));
&G4VEmAdjointModel::DiffCrossSectionFunction2,minEProj+dE1,std::min(minEProj+dE2,maxEProj), 5);
log_ESec_vector->push_back(std::log(std::min(dE2,maxEProj-minEProj)));
log_Prob_vector->push_back(std::log(int_cross_section));
dE1=dE2;
@@ -630,14 +424,11 @@ G4double G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(size_t MatrixIndex,G4
G4AdjointCSMatrix* theMatrix= (*pOnCSMatrixForProdToProjBackwardScattering)[MatrixIndex];
if (IsScatProjToProjCase) theMatrix= (*pOnCSMatrixForScatProjToProjBackwardScattering)[MatrixIndex];
std::vector< G4double >* theLogPrimEnergyVector = theMatrix->GetLogPrimEnergyVector();
//G4double dLog = theMatrix->GetDlog();
std::vector< double>* theLogPrimEnergyVector = theMatrix->GetLogPrimEnergyVector();
if (theLogPrimEnergyVector->size() ==0){
G4cout<<"No data are contained in the given AdjointCSMatrix!"<<std::endl;
G4cout<<"The sampling procedure will be stopped."<<std::endl;
G4cout<<"No data are contained in the given AdjointCSMatrix!"<<G4endl;
G4cout<<"The sampling procedure will be stopped."<<G4endl;
return 0.;
}
@@ -649,11 +440,11 @@ G4double G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(size_t MatrixIndex,G4
G4double aLogPrimEnergy1,aLogPrimEnergy2;
G4double aLogCS1,aLogCS2;
G4double log01,log02;
std::vector< G4double>* aLogSecondEnergyVector1 =0;
std::vector< G4double>* aLogSecondEnergyVector2 =0;
std::vector< G4double>* aLogProbVector1=0;
std::vector< G4double>* aLogProbVector2=0;
G4double log01,log02;
std::vector< double>* aLogSecondEnergyVector1 =0;
std::vector< double>* aLogSecondEnergyVector2 =0;
std::vector< double>* aLogProbVector1=0;
std::vector< double>* aLogProbVector2=0;
std::vector< size_t>* aLogProbVectorIndex1=0;
std::vector< size_t>* aLogProbVectorIndex2=0;
@@ -673,84 +464,31 @@ G4double G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(size_t MatrixIndex,G4
G4double Emin=0.;
G4double Emax=0.;
if (theMatrix->IsScatProjToProjCase()){ //case where Tcut plays a role
//G4cout<<"Here "<<std::endl;
if (ApplyCutInRange) {
Emin=GetSecondAdjEnergyMinForScatProjToProjCase(aPrimEnergy,currentTcutForDirectSecond);
Emax=GetSecondAdjEnergyMaxForScatProjToProjCase(aPrimEnergy);
G4double dE=0;
if (Emin < Emax ){
if (ApplyCutInRange) {
if (second_part_of_same_type && currentTcutForDirectSecond>aPrimEnergy) return aPrimEnergy;
/*if (IsIonisation){
G4double inv_Tcut= 1./currentTcutForDirectSecond;
G4double inv_dE=inv_Tcut-rand_var*(inv_Tcut-1./aPrimEnergy);
Esec= aPrimEnergy+1./inv_dE;
//return Esec;
G4double dE1=currentTcutForDirectSecond;
G4double dE2=currentTcutForDirectSecond*1.00001;
G4double dCS1=DiffCrossSectionMoller(aPrimEnergy+dE1,dE1);
G4double dCS2=DiffCrossSectionMoller(aPrimEnergy+dE2,dE2);
G4double alpha1=std::log(dCS1/dCS2)/std::log(dE1/dE2);
G4double a1=dCS1/std::pow(dE1,alpha1);
dCS1=DiffCrossSectionMoller(aPrimEnergy+dE1,dE1);
dCS2=DiffCrossSectionMoller(aPrimEnergy+dE2,dE2);
return Esec;
dE1=aPrimEnergy/1.00001;
dE2=aPrimEnergy;
dCS1=DiffCrossSectionMoller(aPrimEnergy+dE1,dE1);
dCS2=DiffCrossSectionMoller(aPrimEnergy+dE2,dE2);
G4double alpha2=std::log(dCS1/dCS2)/std::log(dE1/dE2);
G4double a2=dCS1/std::pow(dE1,alpha1);
return Esec;
}*/
log_rand_var1=log_rand_var+theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector1,*aLogProbVector1);
log_rand_var2=log_rand_var+theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector2,*aLogProbVector2);
}
log_dE1 = theInterpolator->Interpolate(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,"Lin");
log_dE2 = theInterpolator->Interpolate(log_rand_var2,*aLogProbVector2,*aLogSecondEnergyVector2,"Lin");
}
log_dE1 = theInterpolator->Interpolate(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,"Lin");
log_dE2 = theInterpolator->Interpolate(log_rand_var2,*aLogProbVector2,*aLogSecondEnergyVector2,"Lin");
dE=std::exp(theInterpolator->LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_dE1,log_dE2));
}
/*log_dE1 = theInterpolator->InterpolateWithIndexVector(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,*aLogProbVectorIndex1,log01,dLog);
log_dE2 = theInterpolator->InterpolateWithIndexVector(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,*aLogProbVectorIndex1,log02,dLog);
*/
Esec = aPrimEnergy +
std::exp(theInterpolator->LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_dE1,log_dE2));
Emin=GetSecondAdjEnergyMinForScatProjToProjCase(aPrimEnergy);
Emax=GetSecondAdjEnergyMaxForScatProjToProjCase(aPrimEnergy);
Esec = aPrimEnergy +dE;
Esec=std::max(Esec,Emin);
Esec=std::min(Esec,Emax);
//G4cout<<"Esec "<<Esec<<std::endl;
//if (Esec > 2.*aPrimEnergy && second_part_of_same_type) Esec = 2.*aPrimEnergy;
}
else { //Tcut condition is already full-filled
/*G4cout<<"Start "<<std::endl;
G4cout<<std::exp((*aLogProbVector1)[0])<<std::endl;
G4cout<<std::exp((*aLogProbVector2)[0])<<std::endl;*/
/*G4double inv_E1= .5/aPrimEnergy;
G4double inv_E=inv_E1-rand_var*(inv_E1-0.00001);
Esec= 1./inv_E;
return Esec;*/
log_E1 = theInterpolator->Interpolate(log_rand_var,*aLogProbVector1,*aLogSecondEnergyVector1,"Lin");
log_E2 = theInterpolator->Interpolate(log_rand_var,*aLogProbVector2,*aLogSecondEnergyVector2,"Lin");
/*log_E1 = theInterpolator->InterpolateWithIndexVector(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,*aLogProbVectorIndex1,log01,dLog);
log_E2 = theInterpolator->InterpolateWithIndexVector(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,*aLogProbVectorIndex1,log02,dLog);
*/
/*G4cout<<std::exp(log_E1)<<std::endl;
G4cout<<std::exp(log_E2)<<std::endl;*/
Esec = std::exp(theInterpolator->LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_E1,log_E2));
Emin=GetSecondAdjEnergyMinForProdToProjCase(aPrimEnergy);
@@ -766,6 +504,39 @@ G4double G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(size_t MatrixIndex,G4
}
//////////////////////////////////////////////////////////////////////////////
//
G4double G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(G4double aPrimEnergy,G4bool IsScatProjToProjCase)
{ SelectCSMatrix(IsScatProjToProjCase);
return SampleAdjSecEnergyFromCSMatrix(indexOfUsedCrossSectionMatrix, aPrimEnergy, IsScatProjToProjCase);
}
//////////////////////////////////////////////////////////////////////////////
//
void G4VEmAdjointModel::SelectCSMatrix(G4bool IsScatProjToProjCase)
{
indexOfUsedCrossSectionMatrix=0;
if (!UseMatrixPerElement) indexOfUsedCrossSectionMatrix = currentMaterialIndex;
else if (!UseOnlyOneMatrixForAllElements) { //Select Material
std::vector<G4double>* CS_Vs_Element = &CS_Vs_ElementForScatProjToProjCase;
lastCS=lastAdjointCSForScatProjToProjCase;
if ( !IsScatProjToProjCase) {
CS_Vs_Element = &CS_Vs_ElementForProdToProjCase;
lastCS=lastAdjointCSForProdToProjCase;
}
G4double rand_var= G4UniformRand();
G4double SumCS=0.;
size_t ind=0;
for (size_t i=0;i<CS_Vs_Element->size();i++){
SumCS+=(*CS_Vs_Element)[i];
if (rand_var<=SumCS/lastCS){
ind=i;
break;
}
}
indexOfUsedCrossSectionMatrix = currentMaterial->GetElement(ind)->GetIndex();
}
}
//////////////////////////////////////////////////////////////////////////////
//
@@ -813,11 +584,40 @@ G4double G4VEmAdjointModel::SampleAdjSecEnergyFromDiffCrossSectionPerAtom(G4doub
}
return E;
}
////////////////////////////////////////////////////////////////////////////////
//
void G4VEmAdjointModel::CorrectPostStepWeight(G4ParticleChange* fParticleChange,
G4double old_weight,
G4double adjointPrimKinEnergy,
G4double projectileKinEnergy,
G4bool IsScatProjToProjCase)
{
G4double new_weight=old_weight;
G4double w_corr =1./CS_biasing_factor;
w_corr*=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
lastCS=lastAdjointCSForScatProjToProjCase;
if ( !IsScatProjToProjCase) lastCS=lastAdjointCSForProdToProjCase;
if (adjointPrimKinEnergy !=preStepEnergy){ //Is that in all cases needed???
G4double post_stepCS=AdjointCrossSection(currentCouple, adjointPrimKinEnergy
,IsScatProjToProjCase );
w_corr*=post_stepCS/lastCS;
}
new_weight*=w_corr;
//G4cout<<"Post step "<<new_weight<<'\t'<<w_corr<<'\t'<<old_weight<<G4endl;
new_weight*=projectileKinEnergy/adjointPrimKinEnergy;//This is needed due to the biasing of diff CS
//by the factor adjointPrimKinEnergy/projectileKinEnergy
fParticleChange->SetParentWeightByProcess(false);
fParticleChange->SetSecondaryWeightByProcess(false);
fParticleChange->ProposeParentWeight(new_weight);
}
//////////////////////////////////////////////////////////////////////////////
//
@@ -829,7 +629,9 @@ G4double G4VEmAdjointModel::GetSecondAdjEnergyMaxForScatProjToProjCase(G4double
//////////////////////////////////////////////////////////////////////////////
//
G4double G4VEmAdjointModel::GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy,G4double Tcut)
{ return PrimAdjEnergy+Tcut;
{ G4double Emin=PrimAdjEnergy;
if (ApplyCutInRange) Emin=PrimAdjEnergy+Tcut;
return Emin;
}
//////////////////////////////////////////////////////////////////////////////
//
@@ -852,16 +654,18 @@ void G4VEmAdjointModel::DefineCurrentMaterial(const G4MaterialCutsCouple* coupl
currentCoupleIndex = couple->GetIndex();
currentMaterialIndex = currentMaterial->GetIndex();
size_t idx=56;
currentTcutForDirectPrim =0.00000000001;
if (theAdjEquivOfDirectPrimPartDef) {
if (theAdjEquivOfDirectPrimPartDef->GetParticleName() == "adj_gamma") idx = 0;
else if (theAdjEquivOfDirectPrimPartDef->GetParticleName() == "adj_e-") idx = 1;
else if (theAdjEquivOfDirectPrimPartDef->GetParticleName() == "adj_e+") idx = 2;
const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
currentTcutForDirectPrim=(*aVec)[currentCoupleIndex];
if (idx <56){
const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
currentTcutForDirectPrim=(*aVec)[currentCoupleIndex];
}
}
currentTcutForDirectSecond =0.00000000001;
if (theAdjEquivOfDirectPrimPartDef == theAdjEquivOfDirectSecondPartDef) {
currentTcutForDirectSecond = currentTcutForDirectPrim;
}
@@ -872,7 +676,36 @@ void G4VEmAdjointModel::DefineCurrentMaterial(const G4MaterialCutsCouple* coupl
else if (theAdjEquivOfDirectSecondPartDef->GetParticleName() == "adj_e+") idx = 2;
const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
currentTcutForDirectSecond=(*aVec)[currentCoupleIndex];
if (idx <56){
const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
currentTcutForDirectPrim=(*aVec)[currentCoupleIndex];
}
}
}
}
}
////////////////////////////////////////////////////////////////////////////////////////////
//
void G4VEmAdjointModel::SetHighEnergyLimit(G4double aVal)
{ HighEnergyLimit=aVal;
if (theDirectEMModel) theDirectEMModel->SetHighEnergyLimit( aVal);
}
////////////////////////////////////////////////////////////////////////////////////////////
//
void G4VEmAdjointModel::SetLowEnergyLimit(G4double aVal)
{
LowEnergyLimit=aVal;
if (theDirectEMModel) theDirectEMModel->SetLowEnergyLimit( aVal);
}
////////////////////////////////////////////////////////////////////////////////////////////
//
void G4VEmAdjointModel::SetAdjointEquivalentOfDirectPrimaryParticleDefinition(G4ParticleDefinition* aPart)
{
theAdjEquivOfDirectPrimPartDef=aPart;
if (theAdjEquivOfDirectPrimPartDef->GetParticleName() =="adj_e-")
theDirectPrimaryPartDef=G4Electron::Electron();
if (theAdjEquivOfDirectPrimPartDef->GetParticleName() =="adj_gamma")
theDirectPrimaryPartDef=G4Gamma::Gamma();
}
source/processes/electromagnetic/adjoint/src/G4eInverseBremsstrahlung.cc
+7 -2
View File
@@ -23,6 +23,9 @@
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
// $Id: G4eInverseBremsstrahlung.cc,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
// GEANT4 tag $Name: geant4-09-03 $
//
#include "G4eInverseBremsstrahlung.hh"
#include "G4VEmAdjointModel.hh"
#include "G4AdjointBremsstrahlungModel.hh"
@@ -30,9 +33,11 @@
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
G4eInverseBremsstrahlung::G4eInverseBremsstrahlung(G4bool whichScatCase,G4String process_name,G4AdjointBremsstrahlungModel* aBremAdjointModel):
G4VAdjointInverseScattering(process_name,whichScatCase)
G4VAdjointReverseReaction(process_name,whichScatCase)
{theAdjointEMModel = aBremAdjointModel;
theAdjointEMModel->SetSecondPartOfSameType(false);
theAdjointEMModel->SetSecondPartOfSameType(false);
if (IsScatProjToProjCase) SetIntegralMode(true);
else SetIntegralMode(false);
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
source/processes/electromagnetic/adjoint/src/G4eInverseCompton.cc
+8 -2
View File
@@ -23,6 +23,9 @@
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
// $Id: G4eInverseCompton.cc,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
// GEANT4 tag $Name: geant4-09-03 $
//
///////////////////////////////////////////////////////
// File name: G4eInverseCompton
//
@@ -38,9 +41,12 @@
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
G4eInverseCompton::G4eInverseCompton(G4bool whichScatCase,G4String process_name,G4AdjointComptonModel* aComptonAdjointModel):
G4VAdjointInverseScattering(process_name,whichScatCase)
G4VAdjointReverseReaction(process_name,whichScatCase)
{theAdjointEMModel = aComptonAdjointModel;
theAdjointEMModel->SetSecondPartOfSameType(false);
theAdjointEMModel->SetSecondPartOfSameType(false);
SetIntegralMode(false);
/*if (IsScatProjToProjCase) SetIntegralMode(false);
else SetIntegralMode(true); */
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
source/processes/electromagnetic/adjoint/src/G4eInverseIonisation.cc
+6 -1
View File
@@ -23,6 +23,9 @@
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
// $Id: G4eInverseIonisation.cc,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
// GEANT4 tag $Name: geant4-09-03 $
//
///////////////////////////////////////////////////////
// File name: G4eInverseIonisation
//
@@ -37,9 +40,11 @@
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
G4eInverseIonisation::G4eInverseIonisation(G4bool whichScatCase,G4String process_name,G4VEmAdjointModel* aEmAdjointModel):
G4VAdjointInverseScattering(process_name,whichScatCase)
G4VAdjointReverseReaction(process_name,whichScatCase)
{theAdjointEMModel = aEmAdjointModel;
theAdjointEMModel->SetSecondPartOfSameType(true);
SetIntegralMode(true);
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
source/processes/electromagnetic/adjoint/src/G4hInverseIonisation.cc
+51
View File
@@ -0,0 +1,51 @@
//
// ********************************************************************
// * License and Disclaimer *
// * *
// * The Geant4 software is copyright of the Copyright Holders of *
// * the Geant4 Collaboration. It is provided under the terms and *
// * conditions of the Geant4 Software License, included in the file *
// * LICENSE and available at http://cern.ch/geant4/license . These *
// * include a list of copyright holders. *
// * *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work make any representation or warranty, express or implied, *
// * regarding this software system or assume any liability for its *
// * use. Please see the license in the file LICENSE and URL above *
// * for the full disclaimer and the limitation of liability. *
// * *
// * This code implementation is the result of the scientific and *
// * technical work of the GEANT4 collaboration. *
// * By using, copying, modifying or distributing the software (or *
// * any work based on the software) you agree to acknowledge its *
// * use in resulting scientific publications, and indicate your *
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
// $Id: G4hInverseIonisation.cc,v 1.2 2009/11/20 10:31:20 ldesorgh Exp $
// GEANT4 tag $Name: geant4-09-03 $
//
///////////////////////////////////////////////////////
// File name: G4hInverseIonisation
//
// Author: Laurent Desorgher
//
// Creation date: 15.02.2009
//
///////////////////////////////////////////////////////
#include "G4hInverseIonisation.hh"
#include "G4VEmAdjointModel.hh"
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
G4hInverseIonisation::G4hInverseIonisation(G4bool whichScatCase,G4String process_name,G4AdjointhIonisationModel* aEmAdjointModel):
G4VAdjointReverseReaction(process_name,whichScatCase)
{theAdjointEMModel = aEmAdjointModel;
theAdjointEMModel->SetSecondPartOfSameType(false);
SetIntegralMode(true);
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
G4hInverseIonisation::~G4hInverseIonisation(){
}