Context Navigation

source: trunk/source/processes/electromagnetic/standard/src/G4PAIModel.cc@ 1330

Last change on this file since 1330 was 1315, checked in by garnier, 15 years ago
update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 32.2 KB

Rev	Line
[819]	1	//
	2	// ********************************************************************
	3	// * License and Disclaimer *
	4	// * *
	5	// * The Geant4 software is copyright of the Copyright Holders of *
	6	// * the Geant4 Collaboration. It is provided under the terms and *
	7	// * conditions of the Geant4 Software License, included in the file *
	8	// * LICENSE and available at http://cern.ch/geant4/license . These *
	9	// * include a list of copyright holders. *
	10	// * *
	11	// * Neither the authors of this software system, nor their employing *
	12	// * institutes,nor the agencies providing financial support for this *
	13	// * work make any representation or warranty, express or implied, *
	14	// * regarding this software system or assume any liability for its *
	15	// * use. Please see the license in the file LICENSE and URL above *
	16	// * for the full disclaimer and the limitation of liability. *
	17	// * *
	18	// * This code implementation is the result of the scientific and *
	19	// * technical work of the GEANT4 collaboration. *
	20	// * By using, copying, modifying or distributing the software (or *
	21	// * any work based on the software) you agree to acknowledge its *
	22	// * use in resulting scientific publications, and indicate your *
	23	// * acceptance of all terms of the Geant4 Software license. *
	24	// ********************************************************************
	25	//
[1315]	26	// $Id: G4PAIModel.cc,v 1.52 2010/06/03 07:28:39 grichine Exp $
	27	// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
[1055]	28	//
	29	// -------------------------------------------------------------------
	30	//
	31	// GEANT4 Class
[819]	32	// File name: G4PAIModel.cc
	33	//
[1315]	34	// Author: Vladimir.Grichine@cern.ch on base of Vladimir Ivanchenko model interface
[819]	35	//
	36	// Creation date: 05.10.2003
	37	//
	38	// Modifications:
	39	//
	40	// 17.08.04 V.Grichine, bug fixed for Tkin<=0 in SampleSecondary
	41	// 16.08.04 V.Grichine, bug fixed in massRatio for DEDX, CrossSection, SampleSecondary
	42	// 08.04.05 Major optimisation of internal interfaces (V.Ivantchenko)
[1196]	43	// 26.07.09 Fixed logic to work with several materials (V.Ivantchenko)
[819]	44	//
	45
	46	#include "G4Region.hh"
	47	#include "G4PhysicsLogVector.hh"
	48	#include "G4PhysicsFreeVector.hh"
	49	#include "G4PhysicsTable.hh"
	50	#include "G4ProductionCutsTable.hh"
	51	#include "G4MaterialCutsCouple.hh"
	52	#include "G4MaterialTable.hh"
	53	#include "G4SandiaTable.hh"
	54	#include "G4OrderedTable.hh"
	55
	56	#include "G4PAIModel.hh"
	57	#include "Randomize.hh"
	58	#include "G4Electron.hh"
	59	#include "G4Positron.hh"
	60	#include "G4Poisson.hh"
	61	#include "G4Step.hh"
	62	#include "G4Material.hh"
	63	#include "G4DynamicParticle.hh"
	64	#include "G4ParticleDefinition.hh"
	65	#include "G4ParticleChangeForLoss.hh"
	66
	67	////////////////////////////////////////////////////////////////////////
	68
	69	using namespace std;
	70
	71	G4PAIModel::G4PAIModel(const G4ParticleDefinition* p, const G4String& nam)
	72	: G4VEmModel(nam),G4VEmFluctuationModel(nam),
	73	fVerbose(0),
	74	fLowestGamma(1.005),
	75	fHighestGamma(10000.),
	76	fTotBin(200),
	77	fMeanNumber(20),
	78	fParticle(0),
	79	fHighKinEnergy(100.*TeV),
	80	fTwoln10(2.0*log(10.0)),
	81	fBg2lim(0.0169),
	82	fTaulim(8.4146e-3)
	83	{
	84	if(p) SetParticle(p);
	85
	86	fElectron = G4Electron::Electron();
	87	fPositron = G4Positron::Positron();
	88
	89	fPAItransferTable = 0;
	90	fPAIdEdxTable = 0;
	91	fSandiaPhotoAbsCof = 0;
	92	fdEdxVector = 0;
	93	fLambdaVector = 0;
	94	fdNdxCutVector = 0;
	95
	96	isInitialised = false;
	97	}
	98
	99	////////////////////////////////////////////////////////////////////////////
	100
	101	G4PAIModel::~G4PAIModel()
	102	{
	103	// G4cout << "PAI: start destruction" << G4endl;
	104	if(fParticleEnergyVector) delete fParticleEnergyVector;
	105	if(fdEdxVector) delete fdEdxVector ;
	106	if(fLambdaVector) delete fLambdaVector;
	107	if(fdNdxCutVector) delete fdNdxCutVector;
	108
	109	if( fPAItransferTable )
	110	{
	111	fPAItransferTable->clearAndDestroy();
	112	delete fPAItransferTable ;
	113	}
	114	if( fPAIdEdxTable )
	115	{
	116	fPAIdEdxTable->clearAndDestroy();
	117	delete fPAIdEdxTable ;
	118	}
	119	if(fSandiaPhotoAbsCof)
	120	{
	121	for(G4int i=0;i<fSandiaIntervalNumber;i++)
	122	{
	123	delete[] fSandiaPhotoAbsCof[i];
	124	}
	125	delete[] fSandiaPhotoAbsCof;
	126	}
	127	//G4cout << "PAI: end destruction" << G4endl;
	128	}
	129
	130	///////////////////////////////////////////////////////////////////////////////
	131
	132	void G4PAIModel::SetParticle(const G4ParticleDefinition* p)
	133	{
	134	if(fParticle == p) return;
	135	fParticle = p;
	136	fMass = fParticle->GetPDGMass();
	137	fSpin = fParticle->GetPDGSpin();
	138	G4double q = fParticle->GetPDGCharge()/eplus;
	139	fChargeSquare = q*q;
	140	fLowKinEnergy = 0.2MeVfMass/proton_mass_c2;
	141	fRatio = electron_mass_c2/fMass;
	142	fQc = fMass/fRatio;
	143	}
	144
	145	////////////////////////////////////////////////////////////////////////////
	146
	147	void G4PAIModel::Initialise(const G4ParticleDefinition* p,
	148	const G4DataVector&)
	149	{
[1315]	150	G4cout<<"G4PAIModel::Initialise for "<<p->GetParticleName()<<G4endl;
[819]	151	if(isInitialised) return;
	152	isInitialised = true;
	153
	154	SetParticle(p);
	155	fLowestKineticEnergy = fMass*(fLowestGamma - 1.0);
	156	fHighestKineticEnergy = fMass*(fHighestGamma - 1.0);
	157
	158	fParticleEnergyVector = new G4PhysicsLogVector(fLowestKineticEnergy,
	159	fHighestKineticEnergy,
	160	fTotBin);
	161
[1055]	162	fParticleChange = GetParticleChangeForLoss();
[819]	163
	164	// Prepare initialization
	165	const G4ProductionCutsTable* theCoupleTable =
	166	G4ProductionCutsTable::GetProductionCutsTable();
	167	const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
	168	size_t numOfMat = G4Material::GetNumberOfMaterials();
	169	size_t numRegions = fPAIRegionVector.size();
	170
	171	for(size_t iReg = 0; iReg < numRegions; ++iReg) // region loop
	172	{
	173	const G4Region* curReg = fPAIRegionVector[iReg];
	174	for(size_t jMat = 0; jMat < numOfMat; ++jMat) // region material loop
	175	{
	176	fMaterial = (*theMaterialTable)[jMat];
	177	fCutCouple = theCoupleTable->GetMaterialCutsCouple( fMaterial,
	178	curReg->GetProductionCuts() );
[1055]	179	//G4cout << "Reg <" <<curReg->GetName() << "> mat <"
	180	// << fMaterial->GetName() << "> fCouple= "
[1196]	181	// << fCutCouple<<" " << p->GetParticleName() <<G4endl;
[819]	182	if( fCutCouple ) {
	183	fMaterialCutsCoupleVector.push_back(fCutCouple);
	184
[1196]	185	fPAItransferTable = new G4PhysicsTable(fTotBin+1);
	186	fPAIdEdxTable = new G4PhysicsTable(fTotBin+1);
	187
[819]	188	fDeltaCutInKinEnergy =
	189	(*theCoupleTable->GetEnergyCutsVector(1))[fCutCouple->GetIndex()];
	190
	191	//ComputeSandiaPhotoAbsCof();
	192	BuildPAIonisationTable();
	193
	194	fPAIxscBank.push_back(fPAItransferTable);
	195	fPAIdEdxBank.push_back(fPAIdEdxTable);
	196	fdEdxTable.push_back(fdEdxVector);
	197
	198	BuildLambdaVector();
	199	fdNdxCutTable.push_back(fdNdxCutVector);
	200	fLambdaTable.push_back(fLambdaVector);
	201	}
	202	}
	203	}
	204	}
	205
	206	//////////////////////////////////////////////////////////////////
	207
[1055]	208	void G4PAIModel::InitialiseMe(const G4ParticleDefinition*)
	209	{}
	210
	211	//////////////////////////////////////////////////////////////////
	212
[819]	213	void G4PAIModel::ComputeSandiaPhotoAbsCof()
	214	{
	215	G4int i, j, numberOfElements ;
	216	static const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
	217
	218	G4SandiaTable thisMaterialSandiaTable(fMatIndex) ;
	219	numberOfElements = (*theMaterialTable)[fMatIndex]->
	220	GetNumberOfElements();
	221	G4int* thisMaterialZ = new G4int[numberOfElements] ;
	222
	223	for(i=0;i<numberOfElements;i++)
	224	{
	225	thisMaterialZ[i] =
	226	(G4int)(*theMaterialTable)[fMatIndex]->GetElement(i)->GetZ() ;
	227	}
	228	fSandiaIntervalNumber = thisMaterialSandiaTable.SandiaIntervals
	229	(thisMaterialZ,numberOfElements) ;
	230
	231	fSandiaIntervalNumber = thisMaterialSandiaTable.SandiaMixing
	232	( thisMaterialZ ,
	233	(*theMaterialTable)[fMatIndex]->GetFractionVector() ,
	234	numberOfElements,fSandiaIntervalNumber) ;
	235
	236	fSandiaPhotoAbsCof = new G4double*[fSandiaIntervalNumber] ;
	237
	238	for(i=0;i<fSandiaIntervalNumber;i++) fSandiaPhotoAbsCof[i] = new G4double[5] ;
	239
	240	for( i = 0 ; i < fSandiaIntervalNumber ; i++ )
	241	{
	242	fSandiaPhotoAbsCof[i][0] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i+1,0) ;
	243
	244	for( j = 1; j < 5 ; j++ )
	245	{
	246	fSandiaPhotoAbsCof[i][j] = thisMaterialSandiaTable.
	247	GetPhotoAbsorpCof(i+1,j)*
	248	(*theMaterialTable)[fMatIndex]->GetDensity() ;
	249	}
	250	}
	251	// delete[] thisMaterialZ ;
	252	}
	253
	254	////////////////////////////////////////////////////////////////////////////
	255	//
	256	// Build tables for the ionization energy loss
	257	// the tables are built for MATERIALS
	258	// *********
	259
	260	void G4PAIModel::BuildPAIonisationTable()
	261	{
	262	G4double LowEdgeEnergy , ionloss ;
	263	G4double tau, Tmax, Tmin, Tkin, deltaLow, gamma, bg2 ;
	264
	265	fdEdxVector = new G4PhysicsLogVector( fLowestKineticEnergy,
	266	fHighestKineticEnergy,
	267	fTotBin);
	268	G4SandiaTable* sandia = fMaterial->GetSandiaTable();
	269
	270	Tmin = sandia->GetSandiaCofForMaterialPAI(0,0)*keV;
	271	deltaLow = 100.eV; // 0.5eV ;
	272
[1196]	273	for (G4int i = 0 ; i <= fTotBin ; i++) //The loop for the kinetic energy
[819]	274	{
	275	LowEdgeEnergy = fParticleEnergyVector->GetLowEdgeEnergy(i) ;
	276	tau = LowEdgeEnergy/fMass ;
	277	gamma = tau +1. ;
	278	// G4cout<<"gamma = "<<gamma<<endl ;
	279	bg2 = tau*( tau + 2. );
	280	Tmax = MaxSecondaryEnergy(fParticle, LowEdgeEnergy);
	281	// Tmax = std::min(fDeltaCutInKinEnergy, Tmax);
	282	Tkin = Tmax ;
	283
	284	if ( Tmax < Tmin + deltaLow ) // low energy safety
	285	Tkin = Tmin + deltaLow ;
	286
	287	fPAIySection.Initialize(fMaterial, Tkin, bg2);
	288
	289	// G4cout<<"ionloss = "<<ionloss*cm/keV<<" keV/cm"<<endl ;
	290	// G4cout<<"n1 = "<<protonPAI.GetIntegralPAIxSection(1)*cm<<" 1/cm"<<endl ;
	291	// G4cout<<"protonPAI.GetSplineSize() = "<<
	292	// protonPAI.GetSplineSize()<<G4endl<<G4endl ;
	293
	294	G4int n = fPAIySection.GetSplineSize();
	295	G4PhysicsFreeVector* transferVector = new G4PhysicsFreeVector(n) ;
	296	G4PhysicsFreeVector* dEdxVector = new G4PhysicsFreeVector(n);
	297
	298	for( G4int k = 0 ; k < n; k++ )
	299	{
	300	transferVector->PutValue( k ,
	301	fPAIySection.GetSplineEnergy(k+1),
	302	fPAIySection.GetIntegralPAIySection(k+1) ) ;
	303	dEdxVector->PutValue( k ,
	304	fPAIySection.GetSplineEnergy(k+1),
	305	fPAIySection.GetIntegralPAIdEdx(k+1) ) ;
	306	}
	307	ionloss = fPAIySection.GetMeanEnergyLoss() ; // total <dE/dx>
	308
	309	if ( ionloss < DBL_MIN) ionloss = DBL_MIN;
	310	fdEdxVector->PutValue(i,ionloss) ;
	311
	312	fPAItransferTable->insertAt(i,transferVector) ;
	313	fPAIdEdxTable->insertAt(i,dEdxVector) ;
	314
	315	} // end of Tkin loop
	316	}
	317
	318	///////////////////////////////////////////////////////////////////////
	319	//
	320	// Build mean free path tables for the delta ray production process
	321	// tables are built for MATERIALS
	322	//
	323
	324	void G4PAIModel::BuildLambdaVector()
	325	{
	326	fLambdaVector = new G4PhysicsLogVector( fLowestKineticEnergy,
	327	fHighestKineticEnergy,
	328	fTotBin ) ;
	329	fdNdxCutVector = new G4PhysicsLogVector( fLowestKineticEnergy,
	330	fHighestKineticEnergy,
	331	fTotBin ) ;
	332	if(fVerbose > 1)
	333	{
	334	G4cout<<"PAIModel DeltaCutInKineticEnergyNow = "
	335	<<fDeltaCutInKinEnergy/keV<<" keV"<<G4endl;
	336	}
[1196]	337	for (G4int i = 0 ; i <= fTotBin ; i++ )
[819]	338	{
	339	G4double dNdxCut = GetdNdxCut(i,fDeltaCutInKinEnergy) ;
[1196]	340	//G4cout << "i= " << i << " x= " << dNdxCut << G4endl;
	341	if(dNdxCut < 0.0) dNdxCut = 0.0;
	342	// G4double lambda = dNdxCut <= DBL_MIN ? DBL_MAX: 1.0/dNdxCut ;
	343	G4double lambda = DBL_MAX;
	344	if(dNdxCut > 0.0) lambda = 1.0/dNdxCut;
[819]	345
	346	fLambdaVector->PutValue(i, lambda) ;
	347	fdNdxCutVector->PutValue(i, dNdxCut) ;
	348	}
	349	}
	350
	351	///////////////////////////////////////////////////////////////////////
	352	//
	353	// Returns integral PAI cross section for energy transfers >= transferCut
	354
	355	G4double
	356	G4PAIModel::GetdNdxCut( G4int iPlace, G4double transferCut)
	357	{
	358	G4int iTransfer;
	359	G4double x1, x2, y1, y2, dNdxCut;
	360	// G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
	361	// G4cout<<"size = "<<G4int((*fPAItransferTable)(iPlace)->GetVectorLength())
	362	// <<G4endl;
	363	for( iTransfer = 0 ;
	364	iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()) ;
	365	iTransfer++)
	366	{
	367	if(transferCut <= (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer))
	368	{
	369	break ;
	370	}
	371	}
	372	if ( iTransfer >= G4int((*fPAItransferTable)(iPlace)->GetVectorLength()) )
	373	{
	374	iTransfer = (*fPAItransferTable)(iPlace)->GetVectorLength() - 1 ;
	375	}
	376	y1 = ((fPAItransferTable)(iPlace))(iTransfer-1) ;
	377	y2 = ((fPAItransferTable)(iPlace))(iTransfer) ;
[1196]	378	if(y1 < 0.0) y1 = 0.0;
	379	if(y2 < 0.0) y2 = 0.0;
[819]	380	// G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
	381	x1 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ;
	382	x2 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
	383	// G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
	384
	385	if ( y1 == y2 ) dNdxCut = y2 ;
	386	else
	387	{
	388	// if ( x1 == x2 ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
[1055]	389	// if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
	390	if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*0.5 ;
[819]	391	else dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ;
	392	}
	393	// G4cout<<""<<dNdxCut<<G4endl;
[1196]	394	if(dNdxCut < 0.0) dNdxCut = 0.0;
[819]	395	return dNdxCut ;
	396	}
	397	///////////////////////////////////////////////////////////////////////
	398	//
	399	// Returns integral dEdx for energy transfers >= transferCut
	400
	401	G4double
	402	G4PAIModel::GetdEdxCut( G4int iPlace, G4double transferCut)
	403	{
	404	G4int iTransfer;
	405	G4double x1, x2, y1, y2, dEdxCut;
[1196]	406	//G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
[819]	407	// G4cout<<"size = "<<G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength())
	408	// <<G4endl;
	409	for( iTransfer = 0 ;
	410	iTransfer < G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength()) ;
	411	iTransfer++)
	412	{
	413	if(transferCut <= (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer))
	414	{
	415	break ;
	416	}
	417	}
	418	if ( iTransfer >= G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength()) )
	419	{
	420	iTransfer = (*fPAIdEdxTable)(iPlace)->GetVectorLength() - 1 ;
	421	}
	422	y1 = ((fPAIdEdxTable)(iPlace))(iTransfer-1) ;
	423	y2 = ((fPAIdEdxTable)(iPlace))(iTransfer) ;
[1196]	424	if(y1 < 0.0) y1 = 0.0;
	425	if(y2 < 0.0) y2 = 0.0;
	426	//G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
[819]	427	x1 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ;
	428	x2 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
[1196]	429	//G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
[819]	430
	431	if ( y1 == y2 ) dEdxCut = y2 ;
	432	else
	433	{
	434	// if ( x1 == x2 ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ;
[1055]	435	// if ( std::abs(x1-x2) <= eV ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ;
	436	if ( std::abs(x1-x2) <= eV ) dEdxCut = y1 + (y2 - y1)*0.5 ;
[819]	437	else dEdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ;
	438	}
[1196]	439	//G4cout<<""<<dEdxCut<<G4endl;
	440	if(dEdxCut < 0.0) dEdxCut = 0.0;
[819]	441	return dEdxCut ;
	442	}
	443
	444	//////////////////////////////////////////////////////////////////////////////
	445
[1055]	446	G4double G4PAIModel::ComputeDEDXPerVolume(const G4Material*,
	447	const G4ParticleDefinition* p,
	448	G4double kineticEnergy,
	449	G4double cutEnergy)
[819]	450	{
	451	G4int iTkin,iPlace;
[1055]	452
	453	//G4double cut = std::min(MaxSecondaryEnergy(p, kineticEnergy), cutEnergy);
	454	G4double cut = cutEnergy;
	455
[819]	456	G4double massRatio = fMass/p->GetPDGMass();
	457	G4double scaledTkin = kineticEnergy*massRatio;
	458	G4double charge = p->GetPDGCharge();
[1055]	459	G4double charge2 = charge*charge;
	460	const G4MaterialCutsCouple* matCC = CurrentCouple();
[819]	461
[1196]	462	size_t jMat = 0;
	463	for(;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
[819]	464	{
	465	if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
	466	}
[1196]	467	//G4cout << "jMat= " << jMat
	468	// << " jMax= " << fMaterialCutsCoupleVector.size()
	469	// << " matCC: " << matCC;
	470	//if(matCC) G4cout << " mat: " << matCC->GetMaterial()->GetName();
	471	// G4cout << G4endl;
	472	if(jMat == fMaterialCutsCoupleVector.size()) return 0.0;
[819]	473
	474	fPAIdEdxTable = fPAIdEdxBank[jMat];
	475	fdEdxVector = fdEdxTable[jMat];
[1196]	476	for(iTkin = 0 ; iTkin <= fTotBin ; iTkin++)
[819]	477	{
	478	if(scaledTkin < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
	479	}
[1196]	480	//G4cout << "E= " << scaledTkin << " iTkin= " << iTkin
	481	// << " " << fdEdxVector->GetVectorLength()<<G4endl;
[819]	482	iPlace = iTkin - 1;
	483	if(iPlace < 0) iPlace = 0;
[1196]	484	else if(iPlace >= fTotBin) iPlace = fTotBin-1;
[1055]	485	G4double dEdx = charge2( (fdEdxVector)(iPlace) - GetdEdxCut(iPlace,cut) );
[819]	486	if( dEdx < 0.) dEdx = 0.;
	487	return dEdx;
	488	}
	489
	490	/////////////////////////////////////////////////////////////////////////
	491
[1055]	492	G4double G4PAIModel::CrossSectionPerVolume( const G4Material*,
	493	const G4ParticleDefinition* p,
	494	G4double kineticEnergy,
	495	G4double cutEnergy,
	496	G4double maxEnergy )
[819]	497	{
	498	G4int iTkin,iPlace;
[1055]	499	G4double tmax = std::min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy);
	500	if(tmax <= cutEnergy) return 0.0;
[819]	501	G4double massRatio = fMass/p->GetPDGMass();
	502	G4double scaledTkin = kineticEnergy*massRatio;
	503	G4double charge = p->GetPDGCharge();
	504	G4double charge2 = charge*charge, cross, cross1, cross2;
[1055]	505	const G4MaterialCutsCouple* matCC = CurrentCouple();
[819]	506
[1196]	507	size_t jMat = 0;
	508	for(;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
[819]	509	{
	510	if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
	511	}
[1196]	512	if(jMat == fMaterialCutsCoupleVector.size()) return 0.0;
[819]	513
	514	fPAItransferTable = fPAIxscBank[jMat];
	515
[1196]	516	for(iTkin = 0 ; iTkin <= fTotBin ; iTkin++)
[819]	517	{
	518	if(scaledTkin < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
	519	}
	520	iPlace = iTkin - 1;
	521	if(iPlace < 0) iPlace = 0;
[1196]	522	else if(iPlace >= fTotBin) iPlace = fTotBin-1;
[819]	523
[1196]	524	//G4cout<<"iPlace = "<<iPlace<<"; tmax = "
[819]	525	// <<tmax<<"; cutEnergy = "<<cutEnergy<<G4endl;
	526	cross1 = GetdNdxCut(iPlace,tmax) ;
[1196]	527	//G4cout<<"cross1 = "<<cross1<<G4endl;
[819]	528	cross2 = GetdNdxCut(iPlace,cutEnergy) ;
[1196]	529	//G4cout<<"cross2 = "<<cross2<<G4endl;
[819]	530	cross = (cross2-cross1)*charge2;
[1196]	531	//G4cout<<"cross = "<<cross<<G4endl;
	532	if( cross < DBL_MIN) cross = 0.0;
[819]	533	// if( cross2 < DBL_MIN) cross2 = DBL_MIN;
	534
	535	// return cross2;
	536	return cross;
	537	}
	538
	539	///////////////////////////////////////////////////////////////////////////
	540	//
	541	// It is analog of PostStepDoIt in terms of secondary electron.
	542	//
	543
	544	void G4PAIModel::SampleSecondaries(std::vector<G4DynamicParticle> vdp,
	545	const G4MaterialCutsCouple* matCC,
	546	const G4DynamicParticle* dp,
	547	G4double tmin,
	548	G4double maxEnergy)
	549	{
[1196]	550	size_t jMat = 0;
	551	for(;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
[819]	552	{
	553	if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
	554	}
[1196]	555	if(jMat == fMaterialCutsCoupleVector.size()) return;
[819]	556
	557	fPAItransferTable = fPAIxscBank[jMat];
	558	fdNdxCutVector = fdNdxCutTable[jMat];
	559
	560	G4double tmax = std::min(MaxSecondaryKinEnergy(dp), maxEnergy);
	561	if( tmin >= tmax && fVerbose > 0)
	562	{
	563	G4cout<<"G4PAIModel::SampleSecondary: tmin >= tmax "<<G4endl;
	564	}
	565	G4ThreeVector direction= dp->GetMomentumDirection();
	566	G4double particleMass = dp->GetMass();
	567	G4double kineticEnergy = dp->GetKineticEnergy();
	568
	569	G4double massRatio = fMass/particleMass;
	570	G4double scaledTkin = kineticEnergy*massRatio;
	571	G4double totalEnergy = kineticEnergy + particleMass;
	572	G4double pSquare = kineticEnergy*(totalEnergy+particleMass);
	573
	574	G4double deltaTkin = GetPostStepTransfer(scaledTkin);
	575
	576	// G4cout<<"G4PAIModel::SampleSecondaries; deltaKIn = "<<deltaTkin/keV<<" keV "<<G4endl;
	577
	578	if( deltaTkin <= 0. && fVerbose > 0)
	579	{
	580	G4cout<<"G4PAIModel::SampleSecondary e- deltaTkin = "<<deltaTkin<<G4endl;
	581	}
	582	if( deltaTkin <= 0.) return;
	583
	584	if( deltaTkin > tmax) deltaTkin = tmax;
	585
	586	G4double deltaTotalMomentum = sqrt(deltaTkin(deltaTkin + 2. electron_mass_c2 ));
	587	G4double totalMomentum = sqrt(pSquare);
	588	G4double costheta = deltaTkin*(totalEnergy + electron_mass_c2)
	589	/(deltaTotalMomentum * totalMomentum);
	590
	591	if( costheta > 0.99999 ) costheta = 0.99999;
	592	G4double sintheta = 0.0;
	593	G4double sin2 = 1. - costheta*costheta;
	594	if( sin2 > 0.) sintheta = sqrt(sin2);
	595
	596	// direction of the delta electron
	597	G4double phi = twopi*G4UniformRand();
	598	G4double dirx = sinthetacos(phi), diry = sinthetasin(phi), dirz = costheta;
	599
	600	G4ThreeVector deltaDirection(dirx,diry,dirz);
	601	deltaDirection.rotateUz(direction);
	602	deltaDirection.unit();
	603
	604	// primary change
	605	kineticEnergy -= deltaTkin;
	606	G4ThreeVector dir = totalMomentumdirection - deltaTotalMomentumdeltaDirection;
	607	direction = dir.unit();
	608	fParticleChange->SetProposedKineticEnergy(kineticEnergy);
	609	fParticleChange->SetProposedMomentumDirection(direction);
	610
	611	// create G4DynamicParticle object for e- delta ray
	612	G4DynamicParticle* deltaRay = new G4DynamicParticle;
	613	deltaRay->SetDefinition(G4Electron::Electron());
	614	deltaRay->SetKineticEnergy( deltaTkin ); // !!! trick for last steps /2.0 ???
	615	deltaRay->SetMomentumDirection(deltaDirection);
	616
	617	vdp->push_back(deltaRay);
	618	}
	619
	620
	621	///////////////////////////////////////////////////////////////////////
	622	//
	623	// Returns post step PAI energy transfer > cut electron energy according to passed
	624	// scaled kinetic energy of particle
	625
	626	G4double
	627	G4PAIModel::GetPostStepTransfer( G4double scaledTkin )
	628	{
	629	// G4cout<<"G4PAIModel::GetPostStepTransfer"<<G4endl ;
	630
	631	G4int iTkin, iTransfer, iPlace ;
	632	G4double transfer = 0.0, position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W ;
	633
[1196]	634	for(iTkin=0;iTkin<=fTotBin;iTkin++)
[819]	635	{
	636	if(scaledTkin < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
	637	}
	638	iPlace = iTkin - 1 ;
	639	// G4cout<<"from search, iPlace = "<<iPlace<<G4endl ;
	640	if(iPlace < 0) iPlace = 0;
[1196]	641	else if(iPlace >= fTotBin) iPlace = fTotBin-1;
[819]	642	dNdxCut1 = (*fdNdxCutVector)(iPlace) ;
	643	// G4cout<<"dNdxCut1 = "<<dNdxCut1<<G4endl ;
	644
	645
	646	if(iTkin == fTotBin) // Fermi plato, try from left
	647	{
	648	position = dNdxCut1*G4UniformRand() ;
	649
	650	for( iTransfer = 0;
	651	iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()); iTransfer++ )
	652	{
	653	if(position >= ((fPAItransferTable)(iPlace))(iTransfer)) break ;
	654	}
	655	transfer = GetEnergyTransfer(iPlace,position,iTransfer);
	656	}
	657	else
	658	{
	659	dNdxCut2 = (*fdNdxCutVector)(iPlace+1) ;
	660	// G4cout<<"dNdxCut2 = "<<dNdxCut2<<G4endl ;
	661	if(iTkin == 0) // Tkin is too small, trying from right only
	662	{
	663	position = dNdxCut2*G4UniformRand() ;
	664
	665	for( iTransfer = 0;
	666	iTransfer < G4int((*fPAItransferTable)(iPlace+1)->GetVectorLength()); iTransfer++ )
	667	{
	668	if(position >= ((fPAItransferTable)(iPlace+1))(iTransfer)) break ;
	669	}
	670	transfer = GetEnergyTransfer(iPlace+1,position,iTransfer);
	671	}
	672	else // general case: Tkin between two vectors of the material
	673	{
	674	E1 = fParticleEnergyVector->GetLowEdgeEnergy(iTkin - 1) ;
	675	E2 = fParticleEnergyVector->GetLowEdgeEnergy(iTkin) ;
	676	W = 1.0/(E2 - E1) ;
	677	W1 = (E2 - scaledTkin)*W ;
	678	W2 = (scaledTkin - E1)*W ;
	679
	680	position = ( dNdxCut1W1 + dNdxCut2W2 )*G4UniformRand() ;
	681
	682	// G4cout<<position<<"\t" ;
	683
	684	G4int iTrMax1, iTrMax2, iTrMax;
	685
	686	iTrMax1 = G4int((*fPAItransferTable)(iPlace)->GetVectorLength());
	687	iTrMax2 = G4int((*fPAItransferTable)(iPlace+1)->GetVectorLength());
	688
	689	if (iTrMax1 >= iTrMax2) iTrMax = iTrMax2;
	690	else iTrMax = iTrMax1;
	691
	692
	693	for( iTransfer = 0; iTransfer < iTrMax; iTransfer++ )
	694	{
	695	if( position >=
	696	( ((fPAItransferTable)(iPlace))(iTransfer)*W1 +
	697	((fPAItransferTable)(iPlace+1))(iTransfer)*W2) ) break ;
	698	}
	699	transfer = GetEnergyTransfer(iPlace,position,iTransfer);
	700	}
	701	}
	702	// G4cout<<"PAImodel PostStepTransfer = "<<transfer/keV<<" keV"<<G4endl ;
	703	if(transfer < 0.0 ) transfer = 0.0 ;
	704	// if(transfer < DBL_MIN ) transfer = DBL_MIN;
	705
	706	return transfer ;
	707	}
	708
	709	///////////////////////////////////////////////////////////////////////
	710	//
	711	// Returns random PAI energy transfer according to passed
	712	// indexes of particle kinetic
	713
	714	G4double
	715	G4PAIModel::GetEnergyTransfer( G4int iPlace, G4double position, G4int iTransfer )
	716	{
	717	G4int iTransferMax;
	718	G4double x1, x2, y1, y2, energyTransfer;
	719
	720	if(iTransfer == 0)
	721	{
	722	energyTransfer = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
	723	}
	724	else
	725	{
	726	iTransferMax = G4int((*fPAItransferTable)(iPlace)->GetVectorLength());
	727
	728	if ( iTransfer >= iTransferMax ) iTransfer = iTransferMax - 1;
	729
	730	y1 = ((fPAItransferTable)(iPlace))(iTransfer-1);
	731	y2 = ((fPAItransferTable)(iPlace))(iTransfer);
	732
	733	x1 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1);
	734	x2 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
	735
	736	if ( x1 == x2 ) energyTransfer = x2;
	737	else
	738	{
	739	if ( y1 == y2 ) energyTransfer = x1 + (x2 - x1)*G4UniformRand();
	740	else
	741	{
	742	energyTransfer = x1 + (position - y1)*(x2 - x1)/(y2 - y1);
	743	}
	744	}
	745	}
	746	return energyTransfer;
	747	}
	748
	749	///////////////////////////////////////////////////////////////////////
	750
	751	G4double G4PAIModel::SampleFluctuations( const G4Material* material,
	752	const G4DynamicParticle* aParticle,
	753	G4double&,
	754	G4double& step,
	755	G4double&)
	756	{
[1196]	757	size_t jMat = 0;
	758	for(;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
[819]	759	{
	760	if( material == fMaterialCutsCoupleVector[jMat]->GetMaterial() ) break;
	761	}
[1196]	762	if(jMat == fMaterialCutsCoupleVector.size()) return 0.0;
[819]	763
	764	fPAItransferTable = fPAIxscBank[jMat];
	765	fdNdxCutVector = fdNdxCutTable[jMat];
	766
[1196]	767	G4int iTkin, iTransfer, iPlace;
[819]	768	G4long numOfCollisions=0;
	769
[1196]	770	//G4cout<<"G4PAIModel::SampleFluctuations"<<G4endl ;
	771	//G4cout<<"in: "<<fMaterialCutsCoupleVector[jMat]->GetMaterial()->GetName()<<G4endl;
[819]	772	G4double loss = 0.0, charge2 ;
	773	G4double stepSum = 0., stepDelta, lambda, omega;
	774	G4double position, E1, E2, W1, W2, W, dNdxCut1, dNdxCut2, meanNumber;
	775	G4bool numb = true;
	776	G4double Tkin = aParticle->GetKineticEnergy() ;
	777	G4double MassRatio = fMass/aParticle->GetDefinition()->GetPDGMass() ;
	778	G4double charge = aParticle->GetDefinition()->GetPDGCharge() ;
	779	charge2 = charge*charge ;
	780	G4double TkinScaled = Tkin*MassRatio ;
	781
[1196]	782	for(iTkin=0;iTkin<=fTotBin;iTkin++)
[819]	783	{
	784	if(TkinScaled < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
	785	}
	786	iPlace = iTkin - 1 ;
	787	if(iPlace < 0) iPlace = 0;
[1196]	788	else if(iPlace >= fTotBin) iPlace = fTotBin - 1;
	789	//G4cout<<"from search, iPlace = "<<iPlace<<G4endl ;
[819]	790	dNdxCut1 = (*fdNdxCutVector)(iPlace) ;
[1196]	791	//G4cout<<"dNdxCut1 = "<<dNdxCut1<<G4endl ;
[819]	792
	793	if(iTkin == fTotBin) // Fermi plato, try from left
	794	{
	795	meanNumber =(((fPAItransferTable)(iPlace))(0)-dNdxCut1)stepcharge2;
	796	if(meanNumber < 0.) meanNumber = 0. ;
	797	// numOfCollisions = RandPoisson::shoot(meanNumber) ;
	798	// numOfCollisions = G4Poisson(meanNumber) ;
	799	if( meanNumber > 0.) lambda = step/meanNumber;
	800	else lambda = DBL_MAX;
	801	while(numb)
	802	{
	803	stepDelta = CLHEP::RandExponential::shoot(lambda);
	804	stepSum += stepDelta;
	805	if(stepSum >= step) break;
	806	numOfCollisions++;
	807	}
[1196]	808	//G4cout<<"##1 numOfCollisions = "<<numOfCollisions<<G4endl ;
[819]	809
	810	while(numOfCollisions)
	811	{
	812	position = dNdxCut1+
	813	(((fPAItransferTable)(iPlace))(0)-dNdxCut1)*G4UniformRand() ;
	814
	815	for( iTransfer = 0;
	816	iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()); iTransfer++ )
	817	{
	818	if(position >= ((fPAItransferTable)(iPlace))(iTransfer)) break ;
	819	}
	820	omega = GetEnergyTransfer(iPlace,position,iTransfer);
[1196]	821	//G4cout<<"G4PAIModel::SampleFluctuations, omega = "<<omega/keV<<" keV; "<<"\t";
[819]	822	loss += omega;
	823	numOfCollisions-- ;
	824	}
	825	}
	826	else
	827	{
	828	dNdxCut2 = (*fdNdxCutVector)(iPlace+1) ;
[1196]	829	//G4cout<<"dNdxCut2 = "<<dNdxCut2<< " iTkin= "<<iTkin<<" iPlace= "<<iPlace<<G4endl;
[819]	830
	831	if(iTkin == 0) // Tkin is too small, trying from right only
	832	{
	833	meanNumber =(((fPAItransferTable)(iPlace+1))(0)-dNdxCut2)stepcharge2;
	834	if( meanNumber < 0. ) meanNumber = 0. ;
	835	// numOfCollisions = CLHEP::RandPoisson::shoot(meanNumber) ;
	836	// numOfCollisions = G4Poisson(meanNumber) ;
[1196]	837	if( meanNumber > 0.) lambda = step/meanNumber;
	838	else lambda = DBL_MAX;
	839	while(numb)
	840	{
	841	stepDelta = CLHEP::RandExponential::shoot(lambda);
	842	stepSum += stepDelta;
	843	if(stepSum >= step) break;
	844	numOfCollisions++;
	845	}
[819]	846
[1196]	847	//G4cout<<"##2 numOfCollisions = "<<numOfCollisions<<G4endl ;
[819]	848
	849	while(numOfCollisions)
	850	{
	851	position = dNdxCut2+
	852	(((fPAItransferTable)(iPlace+1))(0)-dNdxCut2)*G4UniformRand();
	853
	854	for( iTransfer = 0;
	855	iTransfer < G4int((*fPAItransferTable)(iPlace+1)->GetVectorLength()); iTransfer++ )
	856	{
	857	if(position >= ((fPAItransferTable)(iPlace+1))(iTransfer)) break ;
	858	}
	859	omega = GetEnergyTransfer(iPlace,position,iTransfer);
[1196]	860	//G4cout<<omega/keV<<"\t";
[819]	861	loss += omega;
	862	numOfCollisions-- ;
	863	}
	864	}
	865	else // general case: Tkin between two vectors of the material
	866	{
	867	E1 = fParticleEnergyVector->GetLowEdgeEnergy(iTkin - 1) ;
	868	E2 = fParticleEnergyVector->GetLowEdgeEnergy(iTkin) ;
	869	W = 1.0/(E2 - E1) ;
	870	W1 = (E2 - TkinScaled)*W ;
	871	W2 = (TkinScaled - E1)*W ;
	872
[1196]	873	//G4cout << fPAItransferTable->size() << G4endl;
	874	//G4cout<<"((fPAItransferTable)(iPlace))(0) = "<<
[819]	875	// ((fPAItransferTable)(iPlace))(0)<<G4endl ;
[1196]	876	//G4cout<<"((fPAItransferTable)(iPlace+1))(0) = "<<
[819]	877	// ((fPAItransferTable)(iPlace+1))(0)<<G4endl ;
	878
	879	meanNumber=( (((fPAItransferTable)(iPlace))(0)-dNdxCut1)*W1 +
	880	(((fPAItransferTable)(iPlace+1))(0)-dNdxCut2)W2 )step*charge2;
	881	if(meanNumber<0.0) meanNumber = 0.0;
	882	// numOfCollisions = RandPoisson::shoot(meanNumber) ;
	883	// numOfCollisions = G4Poisson(meanNumber) ;
[1196]	884	if( meanNumber > 0.) lambda = step/meanNumber;
	885	else lambda = DBL_MAX;
	886	while(numb)
	887	{
	888	stepDelta = CLHEP::RandExponential::shoot(lambda);
	889	stepSum += stepDelta;
	890	if(stepSum >= step) break;
	891	numOfCollisions++;
	892	}
[819]	893
[1196]	894	//G4cout<<"##3 numOfCollisions = "<<numOfCollisions<<endl ;
[819]	895
	896	while(numOfCollisions)
	897	{
	898	position = dNdxCut1W1 + dNdxCut2W2 +
	899	( ( ((fPAItransferTable)(iPlace))(0)-dNdxCut1 )*W1 +
	900	dNdxCut2+
	901	( ((fPAItransferTable)(iPlace+1))(0)-dNdxCut2 )W2 )G4UniformRand();
	902
	903	// G4cout<<position<<"\t" ;
	904
	905	for( iTransfer = 0;
	906	iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()); iTransfer++ )
	907	{
	908	if( position >=
	909	( ((fPAItransferTable)(iPlace))(iTransfer)*W1 +
	910	((fPAItransferTable)(iPlace+1))(iTransfer)*W2) )
	911	{
	912	break ;
	913	}
	914	}
	915	omega = GetEnergyTransfer(iPlace,position,iTransfer);
	916	// G4cout<<omega/keV<<"\t";
	917	loss += omega;
	918	numOfCollisions-- ;
	919	}
	920	}
	921	}
[1196]	922	//G4cout<<"PAIModel AlongStepLoss = "<<loss/keV<<" keV, on step = "
[819]	923	// <<step/mm<<" mm"<<G4endl ;
	924	if(loss > Tkin) loss=Tkin;
	925	if(loss < 0. ) loss = 0.;
	926	return loss ;
	927
	928	}
	929
	930	//////////////////////////////////////////////////////////////////////
	931	//
	932	// Returns the statistical estimation of the energy loss distribution variance
	933	//
	934
	935
	936	G4double G4PAIModel::Dispersion( const G4Material* material,
	937	const G4DynamicParticle* aParticle,
	938	G4double& tmax,
	939	G4double& step )
	940	{
	941	G4double loss, sumLoss=0., sumLoss2=0., sigma2, meanLoss=0.;
	942	for(G4int i = 0 ; i < fMeanNumber; i++)
	943	{
	944	loss = SampleFluctuations(material,aParticle,tmax,step,meanLoss);
	945	sumLoss += loss;
	946	sumLoss2 += loss*loss;
	947	}
	948	meanLoss = sumLoss/fMeanNumber;
	949	sigma2 = meanLossmeanLoss + (sumLoss2-2sumLoss*meanLoss)/fMeanNumber;
	950	return sigma2;
	951	}
	952
[1055]	953	/////////////////////////////////////////////////////////////////////
[819]	954
[1055]	955	G4double G4PAIModel::MaxSecondaryEnergy( const G4ParticleDefinition* p,
	956	G4double kinEnergy)
	957	{
	958	G4double tmax = kinEnergy;
	959	if(p == fElectron) tmax *= 0.5;
	960	else if(p != fPositron) {
	961	G4double mass = p->GetPDGMass();
	962	G4double ratio= electron_mass_c2/mass;
	963	G4double gamma= kinEnergy/mass + 1.0;
	964	tmax = 2.0electron_mass_c2(gamma*gamma - 1.) /
	965	(1. + 2.0gammaratio + ratio*ratio);
	966	}
	967	return tmax;
	968	}
	969
	970	///////////////////////////////////////////////////////////////
	971
	972	void G4PAIModel::DefineForRegion(const G4Region* r)
	973	{
	974	fPAIRegionVector.push_back(r);
	975	}
	976
[819]	977	//
	978	//
	979	/////////////////////////////////////////////////
	980
	981
	982
	983
	984
	985

Note: See TracBrowser for help on using the repository browser.

Download in other formats: