Context Navigation

source: trunk/source/processes/electromagnetic/lowenergy/src/G4FinalStateIonisationBorn.cc@ 968

Last change on this file since 968 was 961, checked in by garnier, 17 years ago
update processes
Property svn:executable set to ``*
File size: 15.1 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	//
[961]	26	// $Id: G4FinalStateIonisationBorn.cc,v 1.16 2008/12/06 13:47:12 sincerti Exp $
	27	// GEANT4 tag $Name: geant4-09-02-ref-02 $
[819]	28
	29	#include "G4FinalStateIonisationBorn.hh"
	30
[961]	31	//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
[819]	32
	33	G4FinalStateIonisationBorn::G4FinalStateIonisationBorn()
	34	{
	35	G4double scaleFactor = (1.e-22 / 3.343) * m*m;
	36
	37	G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
	38	G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
	39
	40	G4String electron;
	41	G4String proton;
	42
[961]	43	lowEnergyLimitDefault = 12.61 * eV; // SI: i/o 25 eV
[819]	44	highEnergyLimitDefault = 10 * MeV;
	45
	46	char *path = getenv("G4LEDATA");
	47
	48	if (!path)
	49	G4Exception("G4DNACrossSectionDataSet::FullFileName - G4LEDATA environment variable not set");
	50
	51	if (electronDef != 0)
[961]	52	{
	53	electron = electronDef->GetParticleName();
	54	lowEnergyLimit[electron] = 12.61 * eV; // SI: i/o 25 eV
	55	highEnergyLimit[electron] = 30. * keV;
[819]	56
[961]	57	std::ostringstream eFullFileName;
	58	eFullFileName << path << "/dna/sigmadiff_ionisation_e_born.dat";
	59	std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
	60	if (!eDiffCrossSection)
	61	{
	62	G4Exception("G4FinalStateIonisationBorn::ERROR OPENING electron DATA FILE");
	63	}
[819]	64
[961]	65	eTdummyVec.push_back(0.);
	66	while(!eDiffCrossSection.eof())
	67	{
	68	double tDummy;
	69	double eDummy;
	70	eDiffCrossSection>>tDummy>>eDummy;
	71	if (tDummy != eTdummyVec.back()) eTdummyVec.push_back(tDummy);
	72	for (int j=0; j<5; j++)
	73	{
	74	eDiffCrossSection>>eDiffCrossSectionData[j][tDummy][eDummy];
[819]	75
[961]	76	// SI - only if eof is not reached !
	77	if (!eDiffCrossSection.eof()) eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
	78
	79	eVecm[tDummy].push_back(eDummy);
	80
	81	}
[819]	82	}
[961]	83
	84	}
[819]	85	else
[961]	86	{
	87	G4Exception("G4FinalStateIonisationBorn Constructor: electron is not defined");
	88	}
[819]	89
[961]	90	if (protonDef != 0)
	91	{
	92	proton = protonDef->GetParticleName();
	93	lowEnergyLimit[proton] = 500. * keV;
	94	highEnergyLimit[proton] = 10. * MeV;
[819]	95
[961]	96	std::ostringstream pFullFileName;
	97	pFullFileName << path << "/dna/sigmadiff_ionisation_p_born.dat";
	98	std::ifstream pDiffCrossSection(pFullFileName.str().c_str());
	99	if (!pDiffCrossSection)
	100	{
	101	G4Exception("G4FinalStateIonisationBorn::ERROR OPENING proton DATA FILE");
	102	}
	103
	104	pTdummyVec.push_back(0.);
	105	while(!pDiffCrossSection.eof())
[819]	106	{
[961]	107	double tDummy;
	108	double eDummy;
	109	pDiffCrossSection>>tDummy>>eDummy;
	110	if (tDummy != pTdummyVec.back()) pTdummyVec.push_back(tDummy);
	111	for (int j=0; j<5; j++)
	112	{
	113	pDiffCrossSection>>pDiffCrossSectionData[j][tDummy][eDummy];
[819]	114
[961]	115	// SI - only if eof is not reached !
	116	if (!pDiffCrossSection.eof()) pDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
	117
	118	pVecm[tDummy].push_back(eDummy);
	119	}
[819]	120	}
[961]	121	}
[819]	122	else
[961]	123	{
	124	G4Exception("G4FinalStateIonisationBorn Constructor: proton is not defined");
	125	}
[819]	126	}
	127
[961]	128	//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
[819]	129
	130	G4FinalStateIonisationBorn::~G4FinalStateIonisationBorn()
	131	{
	132	eVecm.clear();
	133	pVecm.clear();
	134	}
	135
[961]	136	//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
[819]	137
	138	const G4FinalStateProduct& G4FinalStateIonisationBorn::GenerateFinalState(const G4Track& track, const G4Step& /* step */)
	139	{
	140	product.Clear();
	141
	142	const G4DynamicParticle* particle = track.GetDynamicParticle();
	143
	144	G4double lowLim = lowEnergyLimitDefault;
	145	G4double highLim = highEnergyLimitDefault;
	146
	147	G4double k = particle->GetKineticEnergy();
	148
	149	const G4String& particleName = particle->GetDefinition()->GetParticleName();
	150
	151	std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
	152	pos1 = lowEnergyLimit.find(particleName);
	153
	154	if (pos1 != lowEnergyLimit.end())
[961]	155	{
	156	lowLim = pos1->second;
	157	}
[819]	158
	159	std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
	160	pos2 = highEnergyLimit.find(particleName);
	161
	162	if (pos2 != highEnergyLimit.end())
[961]	163	{
	164	highLim = pos2->second;
	165	}
[819]	166
	167	if (k >= lowLim && k <= highLim)
[961]	168	{
	169	G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
	170	G4double particleMass = particle->GetDefinition()->GetPDGMass();
	171	G4double totalEnergy = k + particleMass;
	172	G4double pSquare = k * (totalEnergy + particleMass);
	173	G4double totalMomentum = std::sqrt(pSquare);
[819]	174
[961]	175	const G4String& particleName = particle->GetDefinition()->GetParticleName();
[819]	176
[961]	177	G4int ionizationShell = cross.RandomSelect(k,particleName);
[819]	178
[961]	179	G4double secondaryKinetic = RandomizeEjectedElectronEnergy(particle->GetDefinition(),k,ionizationShell);
[819]	180
[961]	181	G4double bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
[819]	182
[961]	183	G4double cosTheta = 0.;
	184	G4double phi = 0.;
	185	RandomizeEjectedElectronDirection(track.GetDefinition(), k,secondaryKinetic, cosTheta, phi);
[819]	186
[961]	187	G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
	188	G4double dirX = sinTheta*std::cos(phi);
	189	G4double dirY = sinTheta*std::sin(phi);
	190	G4double dirZ = cosTheta;
	191	G4ThreeVector deltaDirection(dirX,dirY,dirZ);
	192	deltaDirection.rotateUz(primaryDirection);
[819]	193
[961]	194	G4double deltaTotalMomentum = std::sqrt(secondaryKinetic(secondaryKinetic + 2.electron_mass_c2 ));
[819]	195
[961]	196	G4double finalPx = totalMomentumprimaryDirection.x() - deltaTotalMomentumdeltaDirection.x();
	197	G4double finalPy = totalMomentumprimaryDirection.y() - deltaTotalMomentumdeltaDirection.y();
	198	G4double finalPz = totalMomentumprimaryDirection.z() - deltaTotalMomentumdeltaDirection.z();
	199	G4double finalMomentum = std::sqrt(finalPxfinalPx + finalPyfinalPy + finalPz*finalPz);
	200	finalPx /= finalMomentum;
	201	finalPy /= finalMomentum;
	202	finalPz /= finalMomentum;
[819]	203
[961]	204	product.ModifyPrimaryParticle(finalPx,finalPy,finalPz,k-bindingEnergy-secondaryKinetic);
	205	product.AddEnergyDeposit(bindingEnergy);
[819]	206
[961]	207	G4DynamicParticle* aElectron = new G4DynamicParticle(G4Electron::Electron(),deltaDirection,secondaryKinetic);
	208	product.AddSecondary(aElectron);
	209	}
[819]	210
	211	return product;
	212	}
	213
[961]	214	//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
[819]	215
	216	G4double G4FinalStateIonisationBorn::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
[961]	217	G4double k, G4int shell)
[819]	218	{
	219	if (particleDefinition == G4Electron::ElectronDefinition())
[961]	220	{
	221	G4double maximumEnergyTransfer=0.;
	222	if ((k+waterStructure.IonisationEnergy(shell))/2. > k) maximumEnergyTransfer=k;
	223	else maximumEnergyTransfer = (k+waterStructure.IonisationEnergy(shell))/2.;
	224
	225	G4double crossSectionMaximum = 0.;
	226	for(G4double value=waterStructure.IonisationEnergy(shell); value<=maximumEnergyTransfer; value+=0.1*eV)
[819]	227	{
[961]	228	G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
	229	if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
	230	}
[819]	231
[961]	232	G4double secondaryElectronKineticEnergy=0.;
	233	do
	234	{
	235	secondaryElectronKineticEnergy = G4UniformRand() * (maximumEnergyTransfer-waterStructure.IonisationEnergy(shell));
	236	} while(G4UniformRand()*crossSectionMaximum >
	237	DifferentialCrossSection(particleDefinition, k/eV,(secondaryElectronKineticEnergy+waterStructure.IonisationEnergy(shell))/eV,shell));
[819]	238
[961]	239	return secondaryElectronKineticEnergy;
[819]	240
[961]	241	}
[819]	242
	243	if (particleDefinition == G4Proton::ProtonDefinition())
[961]	244	{
	245	G4double maximumKineticEnergyTransfer = 4.* (electron_mass_c2 / proton_mass_c2) * k - (waterStructure.IonisationEnergy(shell));
	246
	247	G4double crossSectionMaximum = 0.;
	248	for (G4double value = waterStructure.IonisationEnergy(shell);
	249	value<=4.*waterStructure.IonisationEnergy(shell) ;
	250	value+=0.1*eV)
[819]	251	{
[961]	252	G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
	253	if (differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
	254	}
[819]	255
[961]	256	G4double secondaryElectronKineticEnergy = 0.;
	257	do
	258	{
	259	secondaryElectronKineticEnergy = G4UniformRand() * maximumKineticEnergyTransfer;
	260	} while(G4UniformRand()*crossSectionMaximum >=
[819]	261	DifferentialCrossSection(particleDefinition, k/eV,(secondaryElectronKineticEnergy+waterStructure.IonisationEnergy(shell))/eV,shell));
	262
[961]	263	return secondaryElectronKineticEnergy;
	264	}
[819]	265
	266	return 0;
	267	}
	268
[961]	269	//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
[819]	270
	271	void G4FinalStateIonisationBorn::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
	272	G4double k,
	273	G4double secKinetic,
	274	G4double & cosTheta,
	275	G4double & phi )
	276	{
	277	if (particleDefinition == G4Electron::ElectronDefinition())
[961]	278	{
	279	phi = twopi * G4UniformRand();
	280	if (secKinetic < 50.eV) cosTheta = (2.G4UniformRand())-1.;
	281	else if (secKinetic <= 200.*eV)
[819]	282	{
[961]	283	if (G4UniformRand() <= 0.1) cosTheta = (2.*G4UniformRand())-1.;
	284	else cosTheta = G4UniformRand()*(std::sqrt(2.)/2);
[819]	285	}
[961]	286	else
	287	{
	288	G4double sin2O = (1.-secKinetic/k) / (1.+secKinetic/(2.*electron_mass_c2));
	289	cosTheta = std::sqrt(1.-sin2O);
	290	}
	291	}
[819]	292
	293	if (particleDefinition == G4Proton::ProtonDefinition())
[961]	294	{
	295	G4double maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
	296	phi = twopi * G4UniformRand();
	297	cosTheta = std::sqrt(secKinetic / maxSecKinetic);
	298	}
[819]	299	}
	300
[961]	301	//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
[819]	302
	303	double G4FinalStateIonisationBorn::DifferentialCrossSection(G4ParticleDefinition * particleDefinition,
	304	G4double k,
	305	G4double energyTransfer,
	306	G4int ionizationLevelIndex)
	307	{
	308	G4double sigma = 0.;
	309
	310	if (energyTransfer >= waterStructure.IonisationEnergy(ionizationLevelIndex))
[961]	311	{
	312	G4double valueT1 = 0;
	313	G4double valueT2 = 0;
	314	G4double valueE21 = 0;
	315	G4double valueE22 = 0;
	316	G4double valueE12 = 0;
	317	G4double valueE11 = 0;
	318
	319	G4double xs11 = 0;
	320	G4double xs12 = 0;
	321	G4double xs21 = 0;
	322	G4double xs22 = 0;
	323
	324	if (particleDefinition == G4Electron::ElectronDefinition())
[819]	325	{
[961]	326	// k should be in eV and energy transfer eV also
[819]	327
[961]	328	std::vector<double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
[819]	329
[961]	330	std::vector<double>::iterator t1 = t2-1;
[819]	331
[961]	332	// SI : the following condition avoids situations where energyTransfer >last vector element
	333	if (energyTransfer <= eVecm[(*t1)].back())
	334	{
	335	std::vector<double>::iterator e12 = std::upper_bound(eVecm[(t1)].begin(),eVecm[(t1)].end(), energyTransfer);
	336	std::vector<double>::iterator e11 = e12-1;
[819]	337
[961]	338	std::vector<double>::iterator e22 = std::upper_bound(eVecm[(t2)].begin(),eVecm[(t2)].end(), energyTransfer);
	339	std::vector<double>::iterator e21 = e22-1;
[819]	340
[961]	341	valueT1 =*t1;
	342	valueT2 =*t2;
	343	valueE21 =*e21;
	344	valueE22 =*e22;
	345	valueE12 =*e12;
	346	valueE11 =*e11;
[819]	347
[961]	348	xs11 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
	349	xs12 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
	350	xs21 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
	351	xs22 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
	352	}
	353
	354	}
[819]	355
[961]	356	if (particleDefinition == G4Proton::ProtonDefinition())
	357	{
	358	// k should be in eV and energy transfer eV also
	359	std::vector<double>::iterator t2 = std::upper_bound(pTdummyVec.begin(),pTdummyVec.end(), k);
	360	std::vector<double>::iterator t1 = t2-1;
	361
	362	std::vector<double>::iterator e12 = std::upper_bound(pVecm[(t1)].begin(),pVecm[(t1)].end(), energyTransfer);
	363	std::vector<double>::iterator e11 = e12-1;
[819]	364
[961]	365	std::vector<double>::iterator e22 = std::upper_bound(pVecm[(t2)].begin(),pVecm[(t2)].end(), energyTransfer);
	366	std::vector<double>::iterator e21 = e22-1;
[819]	367
[961]	368	valueT1 =*t1;
	369	valueT2 =*t2;
	370	valueE21 =*e21;
	371	valueE22 =*e22;
	372	valueE12 =*e12;
	373	valueE11 =*e11;
[819]	374
[961]	375	xs11 = pDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
	376	xs12 = pDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
	377	xs21 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
	378	xs22 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
	379
	380	}
[819]	381
[961]	382	G4double xsProduct = xs11 * xs12 * xs21 * xs22;
	383	if (xsProduct != 0.)
	384	{
	385	sigma = QuadInterpolator( valueE11, valueE12,
[819]	386	valueE21, valueE22,
	387	xs11, xs12,
	388	xs21, xs22,
	389	valueT1, valueT2,
	390	k, energyTransfer);
[961]	391	}
	392
	393	}
	394
[819]	395	return sigma;
	396	}
	397
[961]	398	//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
[819]	399
	400	G4double G4FinalStateIonisationBorn::LogLogInterpolate(G4double e1,
	401	G4double e2,
	402	G4double e,
	403	G4double xs1,
	404	G4double xs2)
	405	{
	406	G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
	407	G4double b = std::log10(xs2) - a*std::log10(e2);
	408	G4double sigma = a*std::log10(e) + b;
	409	G4double value = (std::pow(10.,sigma));
	410	return value;
	411	}
	412
[961]	413	//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
[819]	414
	415	G4double G4FinalStateIonisationBorn::QuadInterpolator(G4double e11, G4double e12,
	416	G4double e21, G4double e22,
	417	G4double xs11, G4double xs12,
	418	G4double xs21, G4double xs22,
	419	G4double t1, G4double t2,
	420	G4double t, G4double e)
	421	{
	422	G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
	423	G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
	424	G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
	425	return value;
	426	}
	427
	428

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

Download in other formats: