Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

source: trunk/source/processes/electromagnetic/muons/test/muEnergyLossTest.cc@ 1228

Last change on this file since 1228 was 1199, checked in by garnier, 16 years ago
nvx fichiers dans CVS
File size: 39.4 KB

Rev	Line
[1199]	1	//
	2	// ********************************************************************
	3	// * License and Disclaimer *
	4	// * *
	5	// * The Geant4 software is copyright of the Copyright Holders of *
	6	// * the Geant4 Collaboration. It is provided under the terms and *
	7	// * conditions of the Geant4 Software License, included in the file *
	8	// * LICENSE and available at http://cern.ch/geant4/license . These *
	9	// * include a list of copyright holders. *
	10	// * *
	11	// * Neither the authors of this software system, nor their employing *
	12	// * institutes,nor the agencies providing financial support for this *
	13	// * work make any representation or warranty, express or implied, *
	14	// * regarding this software system or assume any liability for its *
	15	// * use. Please see the license in the file LICENSE and URL above *
	16	// * for the full disclaimer and the limitation of liability. *
	17	// * *
	18	// * This code implementation is the result of the scientific and *
	19	// * technical work of the GEANT4 collaboration. *
	20	// * By using, copying, modifying or distributing the software (or *
	21	// * any work based on the software) you agree to acknowledge its *
	22	// * use in resulting scientific publications, and indicate your *
	23	// * acceptance of all terms of the Geant4 Software license. *
	24	// ********************************************************************
	25	//
	26	//
	27	// $Id: muEnergyLossTest.cc,v 1.7 2006/06/29 19:49:56 gunter Exp $
	28	// GEANT4 tag $Name: geant4-09-03-cand-01 $
	29	//
	30	//-----------------------------------------------------------------
	31	// new testprogram for testing the mu processes:
	32	// G4MuEnergyLoss
	33	// G4MuIonisation
	34	// G4MuBremstrahlung
	35	// G4MuPairProduction
	36	// G4MuNuclearInteraction
	37	//-----------------------------------------------------------------
	38	// created by L. Urban , May 1998
	39	//---------------------------------------------------------------
	40
	41	#include "G4ios.hh"
	42	#include <fstream>
	43	#include <iomanip>
	44	#include "globals.hh"
	45	#include "G4Timer.hh"
	46	#include "G4MuEnergyLoss.hh"
	47	#include "G4MuIonisation.hh"
	48	#include "G4MuBremsstrahlung.hh"
	49	#include "G4MuPairProduction.hh"
	50	#include "G4MuNuclearInteraction.hh"
	51	#include "G4DynamicParticle.hh"
	52	#include "G4Element.hh"
	53	#include "G4Material.hh"
	54	#include "G4PVPlacement.hh"
	55	#include "G4LogicalVolume.hh"
	56	#include "G4GRSVolume.hh"
	57	#include "G4Box.hh"
	58	#include "G4ProcessManager.hh"
	59	#include "G4Step.hh"
	60	#include "G4StepPoint.hh"
	61	#include "G4Track.hh"
	62	#include "G4Gamma.hh"
	63	#include "G4Electron.hh"
	64	#include "G4Positron.hh"
	65	#include "G4MuonPlus.hh"
	66	#include "G4MuonMinus.hh"
	67	#include "G4GPILSelection.hh"
	68
	69
	70	G4VPhysicalVolume* BuildVolume(G4Material* matworld)
	71	// it builds a simple box filled with material matword .......
	72	{
	73	G4Box myWorldBox= new G4Box ("WBox",10000.cm,10000.cm,10000.cm);
	74
	75	G4LogicalVolume *myWorldLog = new G4LogicalVolume(myWorldBox,matworld,
	76	"WLog",0,0,0) ;
	77
	78	G4PVPlacement *myWorldPhys = new G4PVPlacement(0,G4ThreeVector(),
	79	"WPhys",
	80	myWorldLog,
	81	0,false,0) ;
	82
	83
	84	return myWorldPhys ;
	85
	86	}
	87
	88
	89	int main()
	90	{
	91	//-------- set output format-------
	92	G4cout.setf( std::ios::scientific, std::ios::floatfield );
	93
	94	G4int nrandom;
	95	G4double ran ;
	96
	97	G4cout << "Give the number of random numbers you want to Generate at start!" << G4endl;
	98	G4cin >> nrandom ;
	99	if( nrandom>0)
	100	{
	101	for (G4int ir=0; ir<nrandom; ir++)
	102	{
	103	ran=G4UniformRand() ;
	104	}
	105	}
	106	//--------- Material definition ---------
	107	G4double a, z, ez, density ,temperature,pressure;
	108	G4State state ;
	109	G4String name, symbol;
	110	G4int nel;
	111
	112	a = 9.012*g/mole;
	113	density = 1.848*g/cm3;
	114	G4Material* Be = new G4Material(name="Beryllium", z=4. , a, density);
	115
	116	a = 26.98*g/mole;
	117	density = 2.7*g/cm3;
	118	G4Material* Al = new G4Material(name="Aluminium", z=13., a, density);
	119
	120	a = 28.09*g/mole;
	121	density = 2.33*g/cm3;
	122	G4Material* Si = new G4Material(name="Silicon", z=14., a, density);
	123
	124	G4Element* elH = new G4Element ("Hydrogen", "H", 1. , 1.01*g/mole);
	125
	126	a = 14.01*g/mole;
	127	G4Element* elN = new G4Element(name="Nitrogen", symbol="N", ez=7., a);
	128	a = 16.00*g/mole;
	129	G4Element* elO = new G4Element(name="Oxigen", symbol="O", ez=8., a);
	130	density = 1.29e-03*g/cm3;
	131	state = kStateGas ;
	132	temperature = 273.*kelvin ;
	133	pressure = 1.*atmosphere ;
	134	G4Material* Air = new G4Material(name="Air", density, nel=2 ,
	135	state ,temperature , pressure ) ;
	136	Air->AddElement(elN, .7);
	137	Air->AddElement(elO, .3);
	138
	139	G4Material* H2O = new G4Material ("Water" , 1.*g/cm3, 2);
	140	H2O->AddElement(elH,2);
	141	H2O->AddElement(elO,1);
	142
	143	a = 55.85*g/mole;
	144	density = 7.87*g/cm3;
	145	G4Material* Fe = new G4Material(name="Iron", z=26., a, density);
	146
	147	a = 196.97*g/mole;
	148	density = 19.32*g/cm3;
	149	G4Material* Au = new G4Material(name="Gold", z=79., a, density);
	150
	151	a = 207.19*g/mole;
	152	density = 11.35*g/cm3;
	153	G4Material* Pb = new G4Material(name="Lead", z=82., a, density);
	154
	155	a = 238.03*g/mole;
	156	density = 18.95*g/cm3;
	157	G4Material* U = new G4Material(name="Uranium", z=92., a, density);
	158
	159	//VacuumOnEarth (air with a very low density , like in the beam pipe of an accelerator)
	160	density = 1.0e-10*g/cm3;
	161	G4Material* VacuumOnEarth = new G4Material(name="VacuumOnEarth", density, nel=2);
	162	VacuumOnEarth->AddElement(elN, .7);
	163	VacuumOnEarth->AddElement(elO, .3);
	164
	165	//VacuumInSpace (H , density is extremely small)
	166	a = 1.01*g/mole;
	167	density = 1.0e-50*g/cm3;
	168	G4Material* VacuumInSpace = new G4Material(name="VacuumInSpace",
	169	z=1., a, density) ;
	170
	171	const G4MaterialTable* theMaterialTable ;
	172	G4Material* apttoMaterial ;
	173	G4String MaterialName ;
	174	G4Timer theTimer ;
	175	G4double mloss,sloss,dEdxdelta,dEdxbrems ;
	176	//--------- Particle definition ---------
	177
	178	G4ParticleDefinition* theGamma = G4Gamma::GammaDefinition();
	179	G4ParticleDefinition* theElectron = G4Electron::ElectronDefinition();
	180	G4ParticleDefinition* thePositron = G4Positron::PositronDefinition();
	181	G4ParticleDefinition* theMuonPlus = G4MuonPlus::MuonPlusDefinition();
	182	G4ParticleDefinition* theMuonMinus = G4MuonMinus::MuonMinusDefinition();
	183	G4ParticleDefinition* thePionZero = G4PionZero::PionZeroDefinition();
	184
	185	G4double* GammaKineticEnergyCuts ;
	186	G4double* ElectronKineticEnergyCuts ;
	187	G4double* PositronKineticEnergyCuts ;
	188	G4double* ParticleKineticEnergyCuts ;
	189
	190	theMaterialTable = G4Material::GetMaterialTable() ;
	191
	192	G4double cutinrange,CutInRangeele,CutInRangepos ;
	193
	194	G4ParticleDefinition* theParticle ;
	195	G4GPILSelection selection;
	196
	197	G4double energy, momentum, mass;
	198	G4ProcessVector* palongget ;
	199	G4ProcessVector* palongdo ;
	200	G4ProcessVector* ppostget ;
	201	G4ProcessVector* ppostdo ;
	202
	203	G4String confirm ;
	204
	205	G4cout << " Do you want the mu+ as particle (yes/no)? " << std::flush;
	206	G4cin >> confirm ;
	207	if(confirm == "yes")
	208	{
	209	mass=theMuonPlus->GetPDGMass();
	210	theParticle = theMuonPlus;
	211	}
	212	else
	213	{
	214	G4cout << " Do you want the mu- as particle (yes/no)? " << std::flush;
	215	G4cin >> confirm ;
	216	if(confirm == "yes")
	217	{
	218	mass=theMuonMinus->GetPDGMass();
	219	theParticle = theMuonMinus;
	220	}
	221	}
	222
	223
	224	energy = 1.*GeV + mass ;
	225	momentum=std::sqrt(energyenergy-massmass) ;
	226
	227	G4ParticleMomentum theMomentum(momentum,0.,0.);
	228
	229	G4double pModule = theMomentum.mag();
	230
	231	G4DynamicParticle aParticle(theParticle,energy,theMomentum);
	232
	233	aParticle.SetKineticEnergy(energy-mass);
	234
	235
	236	G4MuIonisation theParticleIonisation ;
	237	G4ProcessManager* theParticleProcessManager = theParticle->GetProcessManager();
	238	theParticleProcessManager->AddProcess(&theParticleIonisation,-1,0,0) ;
	239	G4MuBremsstrahlung theParticleBremsstrahlung ;
	240	theParticleProcessManager->AddProcess(&theParticleBremsstrahlung,-1,-1,1) ;
	241	G4MuPairProduction theParticlePairProduction ;
	242	theParticleProcessManager->AddProcess(&theParticlePairProduction,-1,-1,2) ;
	243	G4MuNuclearInteraction theParticleNuclearInteraction ;
	244	theParticleProcessManager->AddProcess(&theParticleNuclearInteraction,-1,-1,3) ;
	245
	246	G4ForceCondition cond ;
	247	G4ForceCondition* condition = &cond ;
	248
	249	G4double currentSafety ;
	250	G4double& refsafety=currentSafety;
	251
	252	theTimer.Start() ;
	253	G4cout << "cut for GAMMA in mm =" ;
	254	G4cin >> cutinrange ; cutinrange *= mm;
	255	theGamma->SetCuts(cutinrange) ;
	256	G4cout << "gamma,cut in range(mm)=" << theGamma->GetCuts()/mm << G4endl ;
	257
	258	GammaKineticEnergyCuts = theGamma->GetCutsInEnergy() ;
	259	for (G4int icut=0; icut<theMaterialTable->length(); icut++)
	260	{
	261	G4cout << "material index=" << icut << " kin.energy cut(MeV)=" <<
	262	GammaKineticEnergyCuts[icut]/MeV << G4endl ;
	263	}
	264	theTimer.Stop() ;
	265	G4cout << " time = " << theTimer.GetUserElapsed() << G4endl;
	266
	267	theTimer.Start() ;
	268	G4cout << "cut for ELECTRON in mm =" ;
	269	G4cin >> cutinrange ; cutinrange *= mm;
	270	CutInRangeele = cutinrange ;
	271	theElectron->SetCuts(cutinrange) ;
	272	G4cout << "electron,cut in range(mm)=" << theElectron->GetCuts()/mm << G4endl ;
	273
	274	ElectronKineticEnergyCuts = theElectron->GetCutsInEnergy() ;
	275	for ( icut=0; icut<theMaterialTable->length(); icut++)
	276	{
	277	G4cout << "material index=" << icut << " kin.energy cut(MeV)=" <<
	278	ElectronKineticEnergyCuts[icut]/MeV << G4endl ;
	279	}
	280	theTimer.Stop() ;
	281	G4cout << " time = " << theTimer.GetUserElapsed() << G4endl;
	282
	283	theTimer.Start() ;
	284	G4cout << "cut for POSITRON in mm =" ;
	285	G4cin >> cutinrange ; cutinrange *= mm;
	286	CutInRangepos = cutinrange ;
	287	thePositron->SetCuts(cutinrange) ;
	288	G4cout << "positron,cut in range(mm)=" << thePositron->GetCuts()/mm << G4endl ;
	289
	290	PositronKineticEnergyCuts = thePositron->GetCutsInEnergy() ;
	291	for ( icut=0; icut<theMaterialTable->length(); icut++)
	292	{
	293	G4cout << "material index=" << icut << " kin.energy cut(MeV)=" <<
	294	PositronKineticEnergyCuts[icut]/MeV << G4endl ;
	295	}
	296	theTimer.Stop() ;
	297	G4cout << " time = " << theTimer.GetUserElapsed() << G4endl;
	298
	299	theTimer.Start() ;
	300	G4cout << "cut for muons in mm =" ;
	301	G4cin >> cutinrange ; cutinrange *= mm;
	302	theParticle->SetCuts(cutinrange) ;
	303
	304	G4cout << "cut in range(mm)=" << theParticle->GetLengthCuts()/mm << G4endl ;
	305
	306	ParticleKineticEnergyCuts = theParticle->GetEnergyCuts() ;
	307
	308	for ( icut=0; icut<theMaterialTable->length(); icut++)
	309	{
	310	G4cout << "material index=" << icut << " kin.energy cut(MeV)=" <<
	311	ParticleKineticEnergyCuts[icut]/MeV << G4endl ;
	312	}
	313	theTimer.Stop() ;
	314	G4cout << " time = " << theTimer.GetUserElapsed() << G4endl;
	315
	316
	317	G4cout << " ------ ----- " << G4endl ;
	318
	319	palongget = aParticle.GetDefinition()->GetProcessManager()
	320	->GetAlongStepProcessVector(typeGPIL);
	321	ppostget = aParticle.GetDefinition()->GetProcessManager()
	322	->GetPostStepProcessVector(typeGPIL);
	323	palongdo = aParticle.GetDefinition()->GetProcessManager()
	324	->GetAlongStepProcessVector(typeDoIt);
	325	ppostdo = aParticle.GetDefinition()->GetProcessManager()
	326	->GetPostStepProcessVector(typeDoIt);
	327
	328	//---------------------------------- Physics --------------------------------
	329
	330	G4int itry=1, Ntry=1, Nstart, ir;
	331	G4double r ;
	332
	333	//**************************************************************************
	334	const G4int Nbin=137 ;
	335	G4double TkinMeV[Nbin] =
	336	{0.00001,0.000015,0.00002,0.00003,0.00004,0.00005,0.00006,0.00008,
	337	0.0001,0.00015,0.0002,0.0003,0.0004,0.0005,0.0006,0.0008,
	338	0.001,0.0015,0.002,0.003,0.004,0.005,0.006,0.008,
	339	0.01,0.015,0.02,0.03,0.04,0.05,0.06,0.08,
	340	0.1,0.15,0.2,0.3,0.4,0.5,0.6,0.8,
	341	1.,1.5,2.,3.,4.,5.,6.,8.,
	342	10.,15.,20.,30.,40.,50.,60.,80.,
	343	100.,150.,200.,300.,400.,500.,600.,800.,
	344	1.0e3,1.5e3,2.0e3,3.0e3,4.0e3,5.0e3,6.0e3,8.0e3,
	345	1.0e4,1.5e4,2.0e4,3.0e4,4.0e4,5.0e4,6.0e4,8.0e4,
	346	1.0e5,1.5e5,2.0e5,3.0e5,4.0e5,5.0e5,6.0e5,8.0e5,
	347	1.0e6,1.5e6,2.0e6,3.0e6,4.0e6,5.0e6,6.0e6,8.0e6,
	348	1.0e7,1.5e7,2.0e7,3.0e7,4.0e7,5.0e7,6.0e7,8.0e7,
	349	1.0e8,1.5e8,2.0e8,3.0e8,4.0e8,5.0e8,6.0e8,8.0e8,
	350	1.0e9,1.5e9,2.0e9,3.0e9,4.0e9,5.0e9,6.0e9,8.0e9,
	351	1.e10,1.5e10,2.e10,3.e10,4.e10,5.e10,6.e10,8.e10,
	352	1.e11,1.5e11,2.e11,3.e11,4.e11,5.e11,6.e11,8.e11,
	353	1.e12} ;
	354	for (G4int k=0; k<Nbin; k++) TkinMeV[k] *= MeV;
	355
	356	G4int J=-1 ;
	357
	358	G4double lambda,trueStep,geomStep,stepLimit,
	359	previousStepSize,currentMinimumStep ;
	360	G4ParticleChange* aParticleChange ;
	361
	362	G4double T,dEdx,range ;
	363
	364	NEXTMATERIAL: ;
	365	J = J+1 ;
	366	if ( J >= theMaterialTable->length() )
	367	{ G4cout << "that was the last material in the table --> STOP" << G4endl;
	368	return EXIT_FAILURE ; }
	369
	370	apttoMaterial = (*theMaterialTable)[ J ] ;
	371	MaterialName = apttoMaterial->GetName() ;
	372	G4cout << "material=" << MaterialName << G4endl ;
	373	G4cout << "Do you want the Energyloss test 1. for this material?" << G4endl ;
	374	G4cout << "type a positive number if the answer is YES" << G4endl ;
	375	G4cout << "type a negative number if the answer is NO " << G4endl ;
	376	G4int icont ;
	377	G4cin >> icont ;
	378	if ( icont < 0 )
	379	goto NEXTMATERIAL ;
	380
	381	//---------- Volume definition ---------------------
	382
	383	G4VPhysicalVolume* myVolume ;
	384
	385	myVolume = BuildVolume(apttoMaterial) ;
	386
	387	//--------- track and Step definition (for this test ONLY!)------------
	388	G4ThreeVector aPosition(0.,0.,0.);
	389	const G4ThreeVector aDirection(0.,0.,1.) ;
	390	const G4ThreeVector transl(0.,0.,0.) ;
	391
	392	G4double aTime = 0. ;
	393
	394	G4Track* tracke = new G4Track(&aParticle,aTime,aPosition) ;
	395	G4Track& trackele = (*tracke) ;
	396	//(*tracke).SetVolume(myVolume) ;
	397	G4GRSVolume* touche = new G4GRSVolume(myVolume,NULL,transl);
	398	(*tracke).SetTouchable(touche);
	399
	400	(*tracke).SetMomentumDirection(aDirection) ;
	401
	402
	403	G4Step* Step = new G4Step() ;
	404	G4Step& Step = (*Step) ;
	405	tracke->SetStep(Step);
	406
	407	G4StepPoint* aPoint = new G4StepPoint();
	408	(*aPoint).SetPosition(aPosition) ;
	409	G4double safety = 10000.*cm ;
	410	(*aPoint).SetSafety(safety) ;
	411
	412	(*Step).SetPostStepPoint(aPoint) ;
	413
	414	//**************************************************************************
	415
	416	G4cout << G4endl;
	417	G4cout <<" " << MaterialName << " Energyloss test 1." << G4endl ;
	418	G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
	419	G4cout << G4endl ;
	420	G4cout << "kin.en.(MeV) dE/dx(MeV/mm) range(mm) Step(mm)"
	421	<< " dEdx(GeVcm2/g)" << G4endl ;
	422	G4cout << G4endl ;
	423
	424
	425	for ( G4int i=0 ; i<Nbin ; i++)
	426	{
	427	trueStep = cutinrange ;
	428	previousStepSize = cutinrange ;
	429	currentMinimumStep = trueStep ;
	430	(*tracke).SetKineticEnergy(TkinMeV[i]) ;
	431	stepLimit = (*palongget)(0)->AlongStepGetPhysicalInteractionLength(
	432	trackele,
	433	previousStepSize,
	434	currentMinimumStep,
	435	refsafety,
	436	&selection) ;
	437
	438	dEdx = theParticleIonisation.GetdEdx() ;
	439	range = theParticleIonisation.GetRangeNow() ;
	440	T = TkinMeV[i] ;
	441
	442	G4cout <<" " << T/MeV << " " << dEdx/(MeV/mm) << " " ;
	443	G4cout << range/mm << " " << stepLimit/mm << " "
	444	<< dEdx/(GeV/cm)/(apttoMaterial->GetDensity()/(g/cm3)) << G4endl ;
	445
	446	}
	447
	448	G4cout << G4endl;
	449
	450	ENERGYLOSS2: ;
	451
	452	G4cout << "material=" << MaterialName << G4endl ;
	453	G4cout << "Do you want the Energyloss test 2. for this material?" << G4endl ;
	454	G4cout << "type a positive number if the answer is YES" << G4endl ;
	455	G4cout << "type a negative number if the answer is NO " << G4endl ;
	456	G4cin >> icont ;
	457	if ( icont < 0 )
	458	goto ENERGYLOSS3 ;
	459
	460	G4double TMeV,stepmm,stepmx,meanloss,lossnow ;
	461
	462
	463	G4cout << "give an energy value in MeV " ;
	464	G4cin >> TMeV ; TMeV *= MeV;
	465
	466	trueStep = cutinrange ;
	467	previousStepSize = cutinrange ;
	468	currentMinimumStep = trueStep ;
	469	(*tracke).SetKineticEnergy(TMeV) ;
	470	stepmx = (*palongget)(0)->AlongStepGetPhysicalInteractionLength(
	471	trackele,
	472	previousStepSize,
	473	currentMinimumStep,
	474	refsafety,
	475	&selection);
	476
	477	G4cout << " give a steplength in mm , the max. meaningful Step is " << stepmx/mm << " mm" <<G4endl;
	478	G4cout << "Step:" ;
	479	G4cin >> stepmm ; stepmm *= mm;
	480
	481	(*Step).SetTrack(tracke) ;
	482	(*Step).SetStepLength(stepmm);
	483
	484
	485	aParticleChange = (G4ParticleChange*)
	486	((*palongdo)(0)->AlongStepDoIt(trackele,Step));
	487	meanloss = theParticleIonisation.GetMeanLoss() ;
	488	lossnow = TMeV-(*aParticleChange).GetEnergyChange();
	489
	490	G4cout << G4endl;
	491	G4cout <<" " << MaterialName << " Energyloss test 2." << G4endl ;
	492	G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
	493	G4cout << G4endl ;
	494	G4cout << "kin.en.(MeV) Step(mm) meanloss(MeV) act.loss(MeV)" << G4endl ;
	495	G4cout << TMeV/MeV << " " << stepmm/mm << " " << meanloss/MeV << " " << lossnow/MeV << G4endl ;
	496	G4cout << " status change:" << (*aParticleChange).GetStatusChange() << G4endl ;
	497	G4cout << G4endl ;
	498
	499	goto ENERGYLOSS2 ;
	500
	501	ENERGYLOSS3: ;
	502
	503	G4cout << "material=" << MaterialName << G4endl ;
	504	G4cout << "Do you want the Energyloss test 3. for this material?" << G4endl ;
	505	G4cout << "type a positive number if the answer is YES" << G4endl ;
	506	G4cout << "type a negative number if the answer is NO " << G4endl ;
	507	G4cin >> icont ;
	508	if ( icont < 0 )
	509	goto DELTARAY1 ;
	510
	511	G4cout << "give an energy value in MeV " ;
	512	G4cin >> TMeV ; TMeV *= MeV;
	513
	514	trueStep = cutinrange ;
	515	previousStepSize = cutinrange ;
	516	currentMinimumStep = trueStep ;
	517	(*tracke).SetKineticEnergy(TMeV) ;
	518	stepmx = (*palongget)(0)->AlongStepGetPhysicalInteractionLength(
	519	trackele,
	520	previousStepSize,
	521	currentMinimumStep,
	522	refsafety,
	523	&selection);
	524
	525	G4cout << " give a steplength in mm , the max. meaningful Step is " << stepmx << " mm" <<G4endl;
	526	G4cout << "Step:" ;
	527	G4cin >> stepmm ; stepmm *= mm;
	528
	529	(*Step).SetTrack(tracke) ;
	530	(*Step).SetStepLength(stepmm);
	531
	532
	533	G4cout << " give number of events you want " ;
	534	G4int nbev,ibev ;
	535	G4cin >> nbev ;
	536
	537	meanloss=0.;
	538	theTimer.Start();
	539
	540	for ( ibev=0; ibev<nbev; ibev++)
	541	{
	542	aParticleChange = (G4ParticleChange*)
	543	((*palongdo)(0)->AlongStepDoIt(trackele,Step));
	544	lossnow = TMeV-(*aParticleChange).GetEnergyChange();
	545
	546	meanloss += lossnow ;
	547	}
	548
	549	theTimer.Stop();
	550	meanloss /= nbev ;
	551	G4cout << G4endl;
	552	G4cout <<" " << MaterialName << " Energyloss test 3." << G4endl ;
	553	G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
	554	G4cout << G4endl ;
	555	G4cout << "kin.en.(MeV) Step(mm) meanloss(MeV) time/event(sec) " << G4endl ;
	556	G4cout << TMeV/MeV << " " << stepmm/mm << " " << meanloss/MeV << " " <<
	557	theTimer.GetUserElapsed()/nbev << G4endl ;
	558	G4cout << G4endl ;
	559
	560	goto ENERGYLOSS3 ;
	561
	562	DELTARAY1: ;
	563	G4cout << "material=" << MaterialName << G4endl ;
	564	G4cout << "Do you want the delta ray test 1. for this material?" << G4endl ;
	565	G4cout << "type a positive number if the answer is YES" << G4endl ;
	566	G4cout << "type a negative number if the answer is NO " << G4endl ;
	567	G4cin >> icont ;
	568	if ( icont < 0 )
	569	goto DELTARAY2 ;
	570
	571
	572	G4cout << G4endl;
	573	G4cout <<" " << MaterialName << " delta ray test 1." << G4endl ;
	574	G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
	575	G4cout << G4endl ;
	576	G4cout << "kin.en.(MeV) mean free path(mm)" << G4endl ;
	577	G4cout << G4endl ;
	578
	579	for ( i=0 ; i<Nbin ; i++)
	580	{
	581
	582	previousStepSize = cutinrange ;
	583	(*tracke).SetKineticEnergy(TkinMeV[i]) ;
	584	stepLimit = theParticleIonisation.GetMeanFreePath(
	585	trackele,
	586	previousStepSize,
	587	condition) ;
	588
	589	T = TkinMeV[i] ;
	590
	591	G4cout <<" " << T/MeV << " " << stepLimit/mm << G4endl ;
	592
	593	}
	594
	595	G4cout << G4endl;
	596
	597	//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
	598	DELTARAY2: ;
	599
	600	G4cout << "material=" << MaterialName << G4endl ;
	601	G4cout << "Do you want the deltaray test 2. for this material?" << G4endl ;
	602	G4cout << "type a positive number if the answer is YES" << G4endl ;
	603	G4cout << "type a negative number if the answer is NO " << G4endl ;
	604	G4cin >> icont ;
	605
	606	G4double newenergy,dx,dy,dz,Tdelta,ddx,ddy,ddz ;
	607	G4int nd ;
	608	const G4ThreeVector* momdir ;
	609	G4ParticleMomentum ddir ;
	610
	611	if ( icont < 0 )
	612	goto BREMS1 ;
	613
	614	G4cout << "give an energy value in MeV " ;
	615	G4cin >> TMeV ; TMeV *= MeV;
	616
	617	stepmm = 1.*mm ;
	618
	619	(*Step).SetTrack(tracke) ;
	620	(*Step).SetStepLength(stepmm);
	621
	622
	623	(*tracke).SetKineticEnergy(TMeV) ;
	624	aParticleChange = (G4ParticleChange*)
	625	((*ppostdo)(0)->PostStepDoIt(trackele,Step));
	626
	627	newenergy=(*aParticleChange).GetEnergyChange() ;
	628	momdir=(*aParticleChange).GetMomentumChange();
	629	dx = (*momdir).x();
	630	dy = (*momdir).y();
	631	dz = (*momdir).z();
	632	nd=aParticleChange->GetNumberOfSecondaries();
	633
	634	if(nd>0)
	635	{
	636	Tdelta=aParticleChange->GetSecondary(0)->GetKineticEnergy();
	637	ddir=aParticleChange->GetSecondary(0)->
	638	GetMomentumDirection();
	639	ddx = (ddir).x();
	640	ddy = (ddir).y();
	641	ddz = (ddir).z();
	642	}
	643	(*aParticleChange).Clear();
	644
	645
	646	G4cout << G4endl;
	647	G4cout <<" " << MaterialName << " delta ray test 2." << G4endl ;
	648	G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
	649	G4cout << G4endl ;
	650	G4cout << "T=" << TMeV/MeV << " newT=" << newenergy/MeV << " (MeV)" << G4endl ;
	651	G4cout << " status change:" << (*aParticleChange).GetStatusChange() << G4endl ;
	652	if(nd>0)
	653	G4cout << "Tdelta=" << Tdelta/MeV << G4endl ;
	654	G4cout << "new direction:" << dx << " " << dy << " " << dz << G4endl;
	655	if(nd>0)
	656	G4cout << "delta direction:" << ddx << " " << ddy << " " << ddz << G4endl ;
	657	//............................................
	658	nbev=50000 ;
	659	mloss = 0. ;
	660	sloss = 0. ;
	661	theTimer.Start();
	662
	663	for (ibev=0 ; ibev<nbev ; ibev++ )
	664	{
	665	aParticleChange = (G4ParticleChange*)
	666	((*ppostdo)(0)->PostStepDoIt(trackele,Step));
	667	nd=aParticleChange->GetNumberOfSecondaries();
	668	Tdelta=0. ;
	669	if(nd>0)
	670	{
	671	Tdelta=aParticleChange->GetSecondary(0)->GetKineticEnergy();
	672	}
	673	mloss += Tdelta ;
	674	sloss += Tdelta*Tdelta ;
	675	(*aParticleChange).Clear();
	676	}
	677	theTimer.Stop();
	678	mloss /= nbev ;
	679	sloss /= nbev ;
	680	sloss = (sloss-mloss*mloss)/nbev ;
	681	if(sloss>0.)
	682	sloss = std::sqrt(sloss) ;
	683	else
	684	sloss = 0. ;
	685
	686	previousStepSize = cutinrange ;
	687	stepLimit = theParticleIonisation.GetMeanFreePath(
	688	trackele,
	689	previousStepSize,
	690	condition) ;
	691	dEdxdelta=mloss/stepLimit ;
	692
	693	G4cout << " mean energy loss due to delta production (in MeV)=" <<
	694	mloss/MeV << " +- " << sloss/MeV << G4endl ;
	695	G4cout << " dE/dx due to delta production (in MeV/mm,from " <<
	696	nbev << " events )=" << dEdxdelta/(MeV/mm) << G4endl ;
	697	G4cout << "in GeVcm2/g ="
	698	<< dEdxdelta/(GeV/cm)/(apttoMaterial->GetDensity()/(g/cm3)) << G4endl ;
	699	G4cout << " time/delta =" << theTimer.GetUserElapsed()/nbev << G4endl ;
	700	G4cout << G4endl ;
	701
	702	//..................................
	703
	704	goto DELTARAY2 ;
	705
	706	BREMS1: ;
	707	G4cout << "material=" << MaterialName << G4endl ;
	708	G4cout << "Do you want the brems test 1. for this material?" << G4endl ;
	709	G4cout << "type a positive number if the answer is YES" << G4endl ;
	710	G4cout << "type a negative number if the answer is NO " << G4endl ;
	711	G4cin >> icont ;
	712	if ( icont < 0 )
	713	goto BREMS2 ;
	714
	715	G4cout << G4endl;
	716	G4cout <<" " << MaterialName << " brems test 1." << G4endl ;
	717	G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
	718	G4cout << G4endl ;
	719	G4cout << "kin.en.(MeV) mean free path(mm)" << G4endl ;
	720	G4cout << G4endl ;
	721
	722	for ( i=0 ; i<Nbin ; i++)
	723	{
	724
	725	previousStepSize = CutInRangeele ;
	726	(*tracke).SetKineticEnergy(TkinMeV[i]) ;
	727	stepLimit = theParticleBremsstrahlung.GetMeanFreePath(
	728	trackele,
	729	previousStepSize,
	730	condition) ;
	731
	732	T = TkinMeV[i] ;
	733
	734	G4cout <<" " << T/MeV << " " << stepLimit/mm << G4endl ;
	735
	736	}
	737
	738	G4cout << G4endl;
	739
	740	//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
	741	BREMS2: ;
	742
	743	G4cout << "material=" << MaterialName << G4endl ;
	744	G4cout << "Do you want the brems test 2. for this material?" << G4endl ;
	745	G4cout << "type a positive number if the answer is YES" << G4endl ;
	746	G4cout << "type a negative number if the answer is NO " << G4endl ;
	747	G4cin >> icont ;
	748	if ( icont < 0 )
	749	goto PAIR1 ;
	750
	751	G4cout << "give an energy value in MeV " ;
	752	G4cin >> TMeV ; TMeV *= MeV;
	753
	754	stepmm = 1. ;
	755
	756	(*Step).SetTrack(tracke) ;
	757	(*Step).SetStepLength(stepmm);
	758	(*tracke).SetKineticEnergy(TMeV) ;
	759	aParticleChange = (G4ParticleChange*)
	760	((*ppostdo)(1)->PostStepDoIt(trackele,Step));
	761
	762	newenergy=(*aParticleChange).GetEnergyChange() ;
	763	momdir=(*aParticleChange).GetMomentumChange();
	764	dx = (*momdir).x();
	765	dy = (*momdir).y();
	766	dz = (*momdir).z();
	767	nd=aParticleChange->GetNumberOfSecondaries();
	768
	769	if(nd>0)
	770	{
	771	Tdelta=aParticleChange->GetSecondary(0)->GetKineticEnergy();
	772	ddir=aParticleChange->GetSecondary(0)->
	773	GetMomentumDirection();
	774	ddx = (ddir).x();
	775	ddy = (ddir).y();
	776	ddz = (ddir).z();
	777	}
	778	(*aParticleChange).Clear();
	779
	780
	781	G4cout << G4endl;
	782	G4cout <<" " << MaterialName << " brems test 2." << G4endl ;
	783	G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
	784	G4cout << G4endl ;
	785	G4cout << "T=" << TMeV/MeV << " newT=" << newenergy/MeV << " (MeV)" << G4endl ;
	786	G4cout << " status change:" << (*aParticleChange).GetStatusChange() << G4endl ;
	787	if(nd>0)
	788	G4cout << "Tgamma=" << Tdelta/MeV << G4endl ;
	789	G4cout << "new direction:" << dx << " " << dy << " " << dz << G4endl;
	790	if(nd>0)
	791	G4cout << "gamma direction:" << ddx << " " << ddy << " " << ddz << G4endl ;
	792
	793	nbev=50000 ;
	794	mloss = 0. ;
	795	sloss = 0. ;
	796
	797	theTimer.Start();
	798
	799	for (ibev=0 ; ibev<nbev ; ibev++ )
	800	{
	801	aParticleChange = (G4ParticleChange*)
	802	((*ppostdo)(1)->PostStepDoIt(trackele,Step));
	803	nd=aParticleChange->GetNumberOfSecondaries();
	804	Tdelta=0. ;
	805	if(nd>0)
	806	{
	807	Tdelta=aParticleChange->GetSecondary(0)->GetKineticEnergy();
	808	}
	809	mloss += Tdelta ;
	810	sloss += Tdelta*Tdelta ;
	811	(*aParticleChange).Clear();
	812	}
	813	theTimer.Stop();
	814	mloss /= nbev ;
	815	sloss /= nbev ;
	816	sloss = std::sqrt((sloss-mloss*mloss)/nbev) ;
	817
	818	previousStepSize = cutinrange ;
	819	stepLimit = theParticleBremsstrahlung.GetMeanFreePath(
	820	trackele,
	821	previousStepSize,
	822	condition) ;
	823	dEdxbrems=mloss/stepLimit ;
	824
	825	G4cout << " mean energy loss due to bremsstrahlung (in MeV)=" <<
	826	mloss/MeV << " +- " << sloss/MeV << G4endl ;
	827	G4cout << " dE/dx due to bremsstrahlung (in MeV/mm,from " <<
	828	nbev << " events )=" << dEdxbrems/(MeV/mm) << G4endl;
	829	G4cout << "in GeVcm2/g ="
	830	<< dEdxbrems/(GeV/cm)/(apttoMaterial->GetDensity()/(g/cm3)) << G4endl ;
	831	G4cout << " time/brems =" << theTimer.GetUserElapsed()/nbev << G4endl ;
	832	G4cout << G4endl ;
	833
	834	goto BREMS2 ;
	835
	836	PAIR1: ;
	837	G4cout << "material=" << MaterialName << G4endl ;
	838	G4cout << "Do you want the pair test 1. for this material?" << G4endl ;
	839	G4cout << "type a positive number if the answer is YES" << G4endl ;
	840	G4cout << "type a negative number if the answer is NO " << G4endl ;
	841	G4cin >> icont ;
	842	if ( icont < 0 )
	843	goto PAIR2 ;
	844
	845	G4cout << G4endl;
	846	G4cout <<" " << MaterialName << " pair test 1." << G4endl ;
	847	G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
	848	G4cout << G4endl ;
	849	G4cout << "kin.en.(MeV) mean free path(mm)" << G4endl ;
	850	G4cout << G4endl ;
	851
	852	for ( i=0 ; i<Nbin ; i++)
	853	{
	854
	855	previousStepSize = CutInRangeele ;
	856	(*tracke).SetKineticEnergy(TkinMeV[i]) ;
	857	stepLimit = theParticlePairProduction.GetMeanFreePath(
	858	trackele,
	859	previousStepSize,
	860	condition) ;
	861
	862	T = TkinMeV[i] ;
	863
	864	G4cout <<" " << T/MeV << " " << stepLimit/mm << G4endl ;
	865
	866	}
	867
	868	G4cout << G4endl;
	869
	870	//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
	871	PAIR2: ;
	872
	873	G4double Tdelta1,ddx1,ddy1,ddz1,Tdelta2,ddx2,ddy2,ddz2 ;
	874	G4ParticleMomentum ddir1,ddir2 ;
	875	G4double dEdxpair ;
	876
	877	G4cout << "material=" << MaterialName << G4endl ;
	878	G4cout << "Do you want the pair test 2. for this material?" << G4endl ;
	879	G4cout << "type a positive number if the answer is YES" << G4endl ;
	880	G4cout << "type a negative number if the answer is NO " << G4endl ;
	881	G4cin >> icont ;
	882	if ( icont < 0 )
	883	goto NUCL1 ;
	884
	885	G4cout << "give an energy value in MeV " ;
	886	G4cin >> TMeV ; TMeV *= MeV;
	887
	888	stepmm = 1. ;
	889
	890	(*Step).SetTrack(tracke) ;
	891	(*Step).SetStepLength(stepmm);
	892	(*tracke).SetKineticEnergy(TMeV) ;
	893	aParticleChange = (G4ParticleChange*)
	894	((*ppostdo)(2)->PostStepDoIt(trackele,Step));
	895
	896	newenergy=(*aParticleChange).GetEnergyChange() ;
	897	momdir=(*aParticleChange).GetMomentumChange();
	898	dx = (*momdir).x();
	899	dy = (*momdir).y();
	900	dz = (*momdir).z();
	901	nd=aParticleChange->GetNumberOfSecondaries();
	902
	903	if(nd>0)
	904	{
	905	Tdelta1=aParticleChange->GetSecondary(0)->GetKineticEnergy();
	906	ddir1=aParticleChange->GetSecondary(0)->
	907	GetMomentumDirection();
	908	ddx1 = (ddir1).x();
	909	ddy1 = (ddir1).y();
	910	ddz1 = (ddir1).z();
	911	}
	912	if(nd>1)
	913	{
	914	Tdelta2=aParticleChange->GetSecondary(1)->GetKineticEnergy();
	915	ddir2=aParticleChange->GetSecondary(1)->
	916	GetMomentumDirection();
	917	ddx2 = (ddir2).x();
	918	ddy2 = (ddir2).y();
	919	ddz2 = (ddir2).z();
	920	}
	921
	922	(*aParticleChange).Clear();
	923
	924
	925	G4cout << G4endl;
	926	G4cout <<" " << MaterialName << " pair test 2." << G4endl ;
	927	G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
	928	G4cout << G4endl ;
	929	G4cout << "T=" << TMeV/MeV << " newT=" << newenergy/MeV << " (MeV)" << G4endl ;
	930	G4cout << " status change:" << (*aParticleChange).GetStatusChange() << G4endl ;
	931	if(nd>0)
	932	G4cout << "T1=" << Tdelta1/MeV << G4endl ;
	933	if(nd>1)
	934	G4cout << "T2=" << Tdelta2/MeV << G4endl ;
	935
	936	G4cout << "new direction:" << dx << " " << dy << " " << dz << G4endl;
	937	if(nd>0)
	938	G4cout << "direction1:" << ddx1 << " " << ddy1 << " " << ddz1 << G4endl ;
	939	if(nd>1)
	940	G4cout << "direction2:" << ddx2 << " " << ddy2 << " " << ddz2 << G4endl ;
	941
	942	nbev=50000 ;
	943	mloss = 0. ;
	944	sloss = 0. ;
	945
	946	theTimer.Start();
	947
	948	for (ibev=0 ; ibev<nbev ; ibev++ )
	949	{
	950	aParticleChange = (G4ParticleChange*)
	951	((*ppostdo)(2)->PostStepDoIt(trackele,Step));
	952	nd=aParticleChange->GetNumberOfSecondaries();
	953	Tdelta=0. ;
	954	if(nd>0)
	955	{
	956	Tdelta1=aParticleChange->GetSecondary(0)->GetKineticEnergy();
	957	Tdelta=Tdelta1;
	958	}
	959	if(nd>1)
	960	{
	961	Tdelta2=aParticleChange->GetSecondary(1)->GetKineticEnergy();
	962	Tdelta+=Tdelta2;
	963	}
	964
	965	mloss += Tdelta ;
	966	sloss += Tdelta*Tdelta ;
	967	(*aParticleChange).Clear();
	968	}
	969	theTimer.Stop();
	970	mloss /= nbev ;
	971	sloss /= nbev ;
	972	sloss = std::sqrt((sloss-mloss*mloss)/nbev) ;
	973
	974	previousStepSize = cutinrange ;
	975	stepLimit = theParticlePairProduction.GetMeanFreePath(
	976	trackele,
	977	previousStepSize,
	978	condition) ;
	979	dEdxpair=mloss/stepLimit ;
	980
	981	G4cout << " mean energy loss due to pair production (in MeV)=" <<
	982	mloss/MeV << " +- " << sloss/MeV << G4endl ;
	983	G4cout << " dE/dx due to pair production (in MeV/mm,from " <<
	984	nbev << " events )=" << dEdxpair/(MeV/mm) << G4endl;
	985	G4cout << "in GeVcm2/g ="
	986	<< dEdxpair/(GeV/cm)/(apttoMaterial->GetDensity()/(g/cm3)) << G4endl ;
	987	G4cout << " time/pair =" << theTimer.GetUserElapsed()/nbev << G4endl ;
	988	G4cout << G4endl ;
	989
	990	goto PAIR2 ;
	991
	992	NUCL1: ;
	993	G4cout << "material=" << MaterialName << G4endl ;
	994	G4cout << "Do you want the nucl.int. test 1. for this material?" << G4endl ;
	995	G4cout << "type a positive number if the answer is YES" << G4endl ;
	996	G4cout << "type a negative number if the answer is NO " << G4endl ;
	997	G4cin >> icont ;
	998	if ( icont < 0 )
	999	goto NUCL2 ;
	1000
	1001	G4cout << G4endl;
	1002	G4cout <<" " << MaterialName << " nucl.int. test 1." << G4endl ;
	1003	G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
	1004	G4cout << G4endl ;
	1005	G4cout << "kin.en.(MeV) mean free path(mm)" << G4endl ;
	1006	G4cout << G4endl ;
	1007
	1008	for ( i=0 ; i<Nbin ; i++)
	1009	{
	1010
	1011	previousStepSize = CutInRangeele ;
	1012	(*tracke).SetKineticEnergy(TkinMeV[i]) ;
	1013	stepLimit = theParticleNuclearInteraction.GetMeanFreePath(
	1014	trackele,
	1015	previousStepSize,
	1016	condition) ;
	1017
	1018	T = TkinMeV[i] ;
	1019
	1020	G4cout <<" " << T/MeV << " " << stepLimit/mm << G4endl ;
	1021
	1022	}
	1023
	1024	G4cout << G4endl;
	1025
	1026	//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
	1027	NUCL2: ;
	1028
	1029	G4cout << "material=" << MaterialName << G4endl ;
	1030	G4cout << "Do you want the nucl.int. test 2. for this material?" << G4endl ;
	1031	G4cout << "type a positive number if the answer is YES" << G4endl ;
	1032	G4cout << "type a negative number if the answer is NO " << G4endl ;
	1033	G4cin >> icont ;
	1034	if ( icont < 0 )
	1035	goto NEXTMATERIAL ;
	1036
	1037	G4cout << "give an energy value in MeV " ;
	1038	G4cin >> TMeV ; TMeV *= MeV;
	1039
	1040	stepmm = 1. ;
	1041
	1042	(*Step).SetTrack(tracke) ;
	1043	(*Step).SetStepLength(stepmm);
	1044	(*tracke).SetKineticEnergy(TMeV) ;
	1045	aParticleChange = (G4ParticleChange*)
	1046	((*ppostdo)(3)->PostStepDoIt(trackele,Step));
	1047
	1048	newenergy=(*aParticleChange).GetEnergyChange() ;
	1049	momdir=(*aParticleChange).GetMomentumChange();
	1050	dx = (*momdir).x();
	1051	dy = (*momdir).y();
	1052	dz = (*momdir).z();
	1053	nd=aParticleChange->GetNumberOfSecondaries();
	1054
	1055	if(nd>0)
	1056	{
	1057	Tdelta1=aParticleChange->GetSecondary(0)->GetKineticEnergy();
	1058	ddir1=aParticleChange->GetSecondary(0)->
	1059	GetMomentumDirection();
	1060	ddx1 = (ddir1).x();
	1061	ddy1 = (ddir1).y();
	1062	ddz1 = (ddir1).z();
	1063	}
	1064
	1065	(*aParticleChange).Clear();
	1066
	1067
	1068	G4cout << G4endl;
	1069	G4cout <<" " << MaterialName << " nucl.int. test 2." << G4endl ;
	1070	G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
	1071	G4cout << G4endl ;
	1072	G4cout << "T=" << TMeV/MeV << " newT=" << newenergy/MeV << " (MeV)" << G4endl ;
	1073	G4cout << " status change:" << (*aParticleChange).GetStatusChange() << G4endl ;
	1074	if(nd>0)
	1075	G4cout << "T1=" << Tdelta1/MeV << G4endl ;
	1076	if(nd>1)
	1077	G4cout << "T2=" << Tdelta2/MeV << G4endl ;
	1078
	1079	G4cout << "new direction:" << dx << " " << dy << " " << dz << G4endl;
	1080	if(nd>0)
	1081	G4cout << "direction1:" << ddx1 << " " << ddy1 << " " << ddz1 << G4endl ;
	1082	if(nd>1)
	1083	G4cout << "direction2:" << ddx2 << " " << ddy2 << " " << ddz2 << G4endl ;
	1084
	1085	nbev=50000 ;
	1086	mloss = 0. ;
	1087	sloss = 0. ;
	1088
	1089	theTimer.Start();
	1090
	1091	for (ibev=0 ; ibev<nbev ; ibev++ )
	1092	{
	1093	aParticleChange = (G4ParticleChange*)
	1094	((*ppostdo)(3)->PostStepDoIt(trackele,Step));
	1095	nd=aParticleChange->GetNumberOfSecondaries();
	1096	Tdelta=0. ;
	1097	// secondary is not generated ....................................
	1098	Tdelta1 = TMeV - (*aParticleChange).GetEnergyChange() ;
	1099	Tdelta=Tdelta1;
	1100	if(nd>0)
	1101	{
	1102	// secondary is not generated ....................................
	1103	//Tdelta1=aParticleChange->GetSecondary(0)->GetKineticEnergy();
	1104	Tdelta=Tdelta1;
	1105	}
	1106
	1107	mloss += Tdelta ;
	1108	sloss += Tdelta*Tdelta ;
	1109	(*aParticleChange).Clear();
	1110	}
	1111	theTimer.Stop();
	1112	mloss /= nbev ;
	1113	sloss /= nbev ;
	1114	sloss = std::sqrt((sloss-mloss*mloss)/nbev) ;
	1115
	1116	previousStepSize = cutinrange ;
	1117	stepLimit = theParticleNuclearInteraction.GetMeanFreePath(
	1118	trackele,
	1119	previousStepSize,
	1120	condition) ;
	1121	dEdxpair=mloss/stepLimit ;
	1122
	1123	G4cout << " mean energy loss due to nucl.int. production (in MeV)=" <<
	1124	mloss/MeV << " +- " << sloss/MeV << G4endl ;
	1125	G4cout << " dE/dx due to nucl.int. production (in MeV/mm,from " <<
	1126	nbev << " events )=" << dEdxpair/(MeV/mm) << G4endl;
	1127	G4cout << "in GeVcm2/g ="
	1128	<< dEdxpair/(GeV/cm)/(apttoMaterial->GetDensity()/(g/cm3)) << G4endl ;
	1129	G4cout << " time/nucl =" << theTimer.GetUserElapsed()/nbev << G4endl ;
	1130	G4cout << G4endl ;
	1131
	1132	goto NUCL2 ;
	1133	if( J < theMaterialTable->length()-1 )
	1134	goto NEXTMATERIAL ;
	1135
	1136	return EXIT_SUCCESS;
	1137	}

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

Download in other formats: