source: trunk/source/processes/electromagnetic/utils/test/testG4EnergyLossTables.cc@ 1340

Last change on this file since 1340 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: 23.9 KB
Line 
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: testG4EnergyLossTables.cc,v 1.7 2006/06/29 19:55:29 gunter Exp $
28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
29//
30//-------------------------------------------------------------------
31#include "G4ios.hh"
32#include <fstream>
33#include <iomanip>
34#include "globals.hh"
35#include "G4Timer.hh"
36#include "G4MultipleScattering.hh"
37#include "G4DynamicParticle.hh"
38#include "G4Element.hh"
39#include "G4Material.hh"
40#include "G4GRSVolume.hh"
41#include "G4PVPlacement.hh"
42#include "G4LogicalVolume.hh"
43#include "G4Box.hh"
44#include "G4ProcessManager.hh"
45#include "G4Step.hh"
46#include "G4StepPoint.hh"
47#include "G4Track.hh"
48#include "G4Proton.hh"
49#include "G4AntiProton.hh"
50#include "G4PionPlus.hh"
51#include "G4PionMinus.hh"
52#include "G4KaonPlus.hh"
53#include "G4KaonMinus.hh"
54#include "G4MuonPlus.hh"
55#include "G4MuonMinus.hh"
56#include "G4Electron.hh"
57#include "G4Positron.hh"
58#include "G4Gamma.hh"
59#include "G4eEnergyLoss.hh"
60#include "G4eIonisation.hh"
61#include "G4eBremsstrahlung.hh"
62#include "G4MuEnergyLoss.hh"
63#include "G4MuIonisation.hh"
64#include "G4MuBremsstrahlung.hh"
65#include "G4MuPairProduction.hh"
66#include "G4hEnergyLoss.hh"
67#include "G4hIonisation.hh"
68#include "G4GPILSelection.hh"
69
70#include "G4EnergyLossTables.hh"
71
72// test for G4EnergyLossTables, based on tunemsc
73// P. Urban
74
75G4VPhysicalVolume* BuildVolume(G4Material* matworld)
76// it builds a simple box filled with material matword .......
77{
78 G4Box *myWorldBox= new G4Box ("WBox",10000.*cm,10000.*cm,10000.*cm);
79
80 G4LogicalVolume *myWorldLog = new G4LogicalVolume(myWorldBox,matworld,
81 "WLog",0,0,0) ;
82
83 G4PVPlacement *myWorldPhys = new G4PVPlacement(0,G4ThreeVector(),
84 "WPhys",
85 myWorldLog,
86 0,false,0) ;
87
88
89 return myWorldPhys ;
90
91}
92
93
94int main()
95{
96 //-------- set output format-------
97 G4cout.setf( std::ios::scientific, std::ios::floatfield );
98 //---write results to the file msc.out-----
99 std::ofstream outFile("msc.out", std::ios::out ) ;
100 outFile.setf( std::ios::scientific, std::ios::floatfield );
101
102 //--------- Material definition ---------
103 G4Timer theTimer ;
104 G4double a, z, ez, density ,temperature,pressure;
105 G4State state ;
106 G4String name, symbol;
107 G4int nel;
108
109 a = 9.012*g/mole;
110 density = 1.848*g/cm3;
111 G4Material* Be = new G4Material(name="Beryllium", z=4. , a, density);
112
113 a = 12.011*g/mole;
114 density = 2.220*g/cm3;
115 G4Material* C = new G4Material(name="Carbon", z=6. , a, density);
116
117 a = 26.98*g/mole;
118 density = 2.7*g/cm3;
119 G4Material* Al = new G4Material(name="Aluminium", z=13., a, density);
120
121 a = 28.09*g/mole;
122 density = 2.33*g/cm3;
123 G4Material* Si = new G4Material(name="Silicon", z=14., a, density);
124
125 a = 55.85*g/mole;
126 density = 7.87*g/cm3;
127 G4Material* Fe = new G4Material(name="Iron", z=26., a, density);
128
129 a = 63.540*g/mole;
130 density = 8.960*g/cm3;
131 G4Material* Cu = new G4Material(name="Copper", z=29., a, density);
132
133 a = 196.97*g/mole;
134 density = 19.32*g/cm3;
135 G4Material* Au = new G4Material(name="Gold", z=79., a, density);
136
137 a = 207.19*g/mole;
138 density = 11.35*g/cm3;
139 G4Material* Pb = new G4Material(name="Lead", z=82., a, density);
140
141 a = 238.03*g/mole;
142 density = 18.7*g/cm3;
143 G4Material* U = new G4Material(name="Uranium", z=92., a, density);
144
145 a = 14.01*g/mole;
146 G4Element* elN = new G4Element(name="Nitrogen", symbol="N", ez=7., a);
147 a = 16.00*g/mole;
148 G4Element* elO = new G4Element(name="Oxigen", symbol="O", ez=8., a);
149 density = 1.29e-03*g/cm3;
150 state = kStateGas ;
151 temperature = 273.*kelvin ;
152 pressure = 1.*atmosphere ;
153 G4Material* Air = new G4Material(name="Air", density, nel=2 ,
154 state ,temperature , pressure ) ;
155 Air->AddElement(elN, .7);
156 Air->AddElement(elO, .3);
157
158 a = 0. ;
159 density = 0. ;
160 G4Material* Vac= new G4Material(name="Vacuum",z=0., a, density,kVacuum);
161
162 static const G4MaterialTable* theMaterialTable ;
163 G4Material* apttoMaterial ;
164 G4String MaterialName ;
165
166//--------- Particle definition ---------
167
168 G4ParticleDefinition* theGamma = G4Gamma::GammaDefinition();
169
170 G4ParticleDefinition* theElectron = G4Electron::ElectronDefinition();
171 G4ParticleDefinition* thePositron = G4Positron::PositronDefinition();
172
173 G4ParticleDefinition* theMuonPlus = G4MuonPlus::MuonPlusDefinition();
174 G4ParticleDefinition* theMuonMinus = G4MuonMinus::MuonMinusDefinition();
175
176 G4ParticleDefinition* theProton = G4Proton::ProtonDefinition();
177 G4ParticleDefinition* theAntiProton = G4AntiProton::AntiProtonDefinition();
178 G4ParticleDefinition* thePionPlus = G4PionPlus::PionPlusDefinition();
179 G4ParticleDefinition* thePionMinus = G4PionMinus::PionMinusDefinition();
180 G4ParticleDefinition* theKaonPlus = G4KaonPlus::KaonPlusDefinition();
181 G4ParticleDefinition* theKaonMinus = G4KaonMinus::KaonMinusDefinition();
182
183//--------- Process definition ---------
184 G4eIonisation theElectronIonisation,thePositronIonisation;
185 G4eBremsstrahlung theElectronBremsstrahlung,thePositronBremsstrahlung;
186//..........e-/e+....................................................
187 G4MultipleScattering theElectronMultipleScattering,thePositronMultipleScattering ;
188
189 G4ProcessManager* theElectronProcessManager = theElectron->GetProcessManager();
190 G4ProcessManager* thePositronProcessManager = thePositron->GetProcessManager();
191
192 theElectronProcessManager->AddProcess(&theElectronMultipleScattering,-1,0,0) ;
193 thePositronProcessManager->AddProcess(&thePositronMultipleScattering,-1,0,0) ;
194 theElectronProcessManager->AddProcess(&theElectronIonisation,-1,1,1) ;
195 theElectronProcessManager->AddProcess(&theElectronBremsstrahlung,-1,-1,2) ;
196 thePositronProcessManager->AddProcess(&thePositronIonisation,-1,1,1) ;
197 thePositronProcessManager->AddProcess(&thePositronBremsstrahlung,-1,-1,2) ;
198
199//..........mu+/mu-................................................
200 G4MuIonisation theMuonPlusIonisation,theMuonMinusIonisation ;
201 G4MuBremsstrahlung theMuonPlusBremsstrahlung,theMuonMinusBremsstrahlung ;
202 G4MuPairProduction theMuonPlusPairProduction,theMuonMinusPairProduction ;
203 G4MultipleScattering theMuonPlusMultipleScattering,theMuonMinusMultipleScattering;
204
205 G4ProcessManager* theMuonPlusProcessManager = theMuonPlus->GetProcessManager();
206 theMuonPlusProcessManager->AddProcess(&theMuonPlusMultipleScattering,-1,0,0) ;
207 theMuonPlusProcessManager->AddProcess(&theMuonPlusIonisation,-1,1,1);
208 theMuonPlusProcessManager->AddProcess(&theMuonPlusBremsstrahlung,-1,-1,2);
209 theMuonPlusProcessManager->AddProcess(&theMuonPlusPairProduction,-1,-1,3);
210
211 G4ProcessManager* theMuonMinusProcessManager = theMuonMinus->GetProcessManager();
212 theMuonMinusProcessManager->AddProcess(&theMuonMinusMultipleScattering,-1,0,0) ;
213 theMuonMinusProcessManager->AddProcess(&theMuonMinusIonisation,-1,1,1);
214 theMuonMinusProcessManager->AddProcess(&theMuonMinusBremsstrahlung,-1,-1,2);
215 theMuonMinusProcessManager->AddProcess(&theMuonMinusPairProduction,-1,-1,3);
216
217//------multiple instantiation of G4MultipleScattering for hadrons-----------------
218 G4MultipleScattering theProtonMultipleScattering,theAntiProtonMultipleScattering;
219 G4hIonisation theProtonIonisation,theAntiProtonIonisation ;
220 G4MultipleScattering thePionPlusMultipleScattering,thePionMinusMultipleScattering;
221 G4hIonisation thePionPlusIonisation,thePionMinusIonisation;
222 G4MultipleScattering theKaonPlusMultipleScattering,theKaonMinusMultipleScattering;
223 G4hIonisation theKaonPlusIonisation,theKaonMinusIonisation;
224
225 G4ProcessManager* theProtonProcessManager = theProton->GetProcessManager();
226 theProtonProcessManager->AddProcess(&theProtonMultipleScattering,-1,0,0) ;
227 theProtonProcessManager->AddProcess(&theProtonIonisation,-1,1,1);
228
229 G4ProcessManager* thePionPlusProcessManager = thePionPlus->GetProcessManager();
230 thePionPlusProcessManager->AddProcess(&thePionPlusMultipleScattering,-1,0,0) ;
231 thePionPlusProcessManager->AddProcess(&thePionPlusIonisation,-1,1,1);
232
233 G4ProcessManager* theKaonPlusProcessManager = theKaonPlus->GetProcessManager();
234 theKaonPlusProcessManager->AddProcess(&theKaonPlusMultipleScattering,-1,0,0) ;
235 theKaonPlusProcessManager->AddProcess(&theKaonPlusIonisation,-1,1,1);
236
237 G4ProcessManager* theAntiProtonProcessManager = theAntiProton->GetProcessManager();
238 theAntiProtonProcessManager->AddProcess(&theAntiProtonMultipleScattering,-1,0,0) ;
239 theAntiProtonProcessManager->AddProcess(&theAntiProtonIonisation,-1,1,1);
240
241 G4ProcessManager* thePionMinusProcessManager = thePionMinus->GetProcessManager();
242 thePionMinusProcessManager->AddProcess(&thePionMinusMultipleScattering,-1,0,0) ;
243 thePionMinusProcessManager->AddProcess(&thePionMinusIonisation,-1,1,1);
244
245 G4ProcessManager* theKaonMinusProcessManager = theKaonMinus->GetProcessManager();
246 theKaonMinusProcessManager->AddProcess(&theKaonMinusMultipleScattering,-1,0,0) ;
247 theKaonMinusProcessManager->AddProcess(&theKaonMinusIonisation,-1,1,1);
248
249 G4GPILSelection selection;
250 G4ParticleWithCuts* theParticle ;
251 G4MultipleScattering* theParticleMultipleScattering;
252 G4ProcessManager* theParticleProcessManager ;
253 G4String confirm ;
254 G4int i1e ;
255 G4double theta1e,th1,th2 ;
256 G4double theta1edata,errth,lambdadata,lambdadatamin,lambdadatamax ;
257
258 // loop until no particle is selected
259 while (1) {
260
261 G4cout << "Do you want the electron as particle (yes/no)?" << std::flush;
262 G4cin >> confirm ;
263 if(confirm == "yes")
264 {
265 theParticle = theElectron ;
266 theParticleMultipleScattering=&theElectronMultipleScattering;
267 theParticleProcessManager=theElectronProcessManager;
268 outFile << " ----------particle = electron -------------" << G4endl;
269 }
270 else
271 {
272 G4cout << "Do you want the positron as particle (yes/no)?" << std::flush;
273 G4cin >> confirm ;
274 if(confirm == "yes")
275 {
276 theParticle = thePositron ;
277 theParticleMultipleScattering=&thePositronMultipleScattering;
278 theParticleProcessManager=thePositronProcessManager;
279 outFile << " ----------particle = positron -------------" << G4endl;
280 }
281 else
282 {
283 G4cout << "Do you want the mu+ as particle (yes/no)?" << std::flush;
284 G4cin >> confirm ;
285 if(confirm == "yes")
286 {
287 theParticle = theMuonPlus ;
288 theParticleMultipleScattering=&theMuonPlusMultipleScattering;
289 theParticleProcessManager=theMuonPlusProcessManager;
290 outFile << " --------particle = mu+ -------------" << G4endl;
291 }
292 else
293 {
294 G4cout << "Do you want the mu- as particle (yes/no)?" << std::flush;
295 G4cin >> confirm ;
296 if(confirm == "yes")
297 {
298 theParticle = theMuonMinus ;
299 theParticleMultipleScattering=&theMuonMinusMultipleScattering;
300 theParticleProcessManager=theMuonMinusProcessManager;
301 outFile << " --------particle = mu- -------------" << G4endl;
302 }
303 else
304 {
305 G4cout << " Do you want the proton as particle (yes/no)? " << std::flush;
306 G4cin >> confirm ;
307 if(confirm == "yes")
308 {
309 theParticle = theProton;
310 theParticleMultipleScattering=&theProtonMultipleScattering;
311 theParticleProcessManager=theProtonProcessManager;
312 outFile << " ---------- particle = proton ----------------" << G4endl;
313 }
314 else
315 {
316 G4cout << " Do you want the antiproton as particle (yes/no)? " << std::flush;
317 G4cin >> confirm ;
318 if(confirm == "yes")
319 {
320 theParticle = theAntiProton;
321 theParticleMultipleScattering=&theAntiProtonMultipleScattering;
322 theParticleProcessManager=theAntiProtonProcessManager;
323 outFile << " ---------- particle = antiproton ----------------" << G4endl;
324 }
325 else
326 {
327 G4cout << " Do you want the pi+ as particle (yes/no)? " << std::flush;
328 G4cin >> confirm ;
329 if(confirm == "yes")
330 {
331 theParticle = thePionPlus;
332 theParticleMultipleScattering=&thePionPlusMultipleScattering;
333 theParticleProcessManager=thePionPlusProcessManager;
334 outFile << " ---------- particle = pi+ ----------------" << G4endl;
335 }
336 else
337 {
338 G4cout << " Do you want the pi- as particle (yes/no)? " << std::flush;
339 G4cin >> confirm ;
340 if(confirm == "yes")
341 {
342 theParticle = thePionMinus;
343 theParticleMultipleScattering=&thePionMinusMultipleScattering;
344 theParticleProcessManager=thePionMinusProcessManager;
345 outFile << " ---------- particle = pi- ----------------" << G4endl;
346 }
347 else
348 {
349 G4cout << " Do you want the K+ as particle (yes/no)? " << std::flush;
350 G4cin >> confirm ;
351 if(confirm == "yes")
352 {
353 theParticle = theKaonPlus;
354 theParticleMultipleScattering=&theKaonPlusMultipleScattering;
355 theParticleProcessManager=theKaonPlusProcessManager;
356 outFile << " ---------- particle = K+ ----------------" << G4endl;
357 }
358 else
359 {
360 G4cout << " Do you want the K- as particle (yes/no)? " << std::flush;
361 G4cin >> confirm ;
362 if(confirm == "yes")
363 {
364 theParticle = theKaonMinus;
365 theParticleMultipleScattering=&theKaonMinusMultipleScattering;
366 theParticleProcessManager=theKaonMinusProcessManager;
367 outFile << " ---------- particle = K- ----------------" << G4endl;
368 }
369 else
370 {
371 G4cout << " There is no other particle in the test." << G4endl;
372 return EXIT_SUCCESS;
373 }
374 }
375 }
376 }
377 }
378 }
379 }
380 }
381 }
382 }
383
384 G4ForceCondition cond;
385 G4ForceCondition* condition=&cond ;
386
387 G4double* ElectronKineticEnergyCuts ;
388 G4double* PositronKineticEnergyCuts ;
389 G4double* GammaKineticEnergyCuts ;
390 G4double* ParticleKineticEnergyCuts ;
391
392 theMaterialTable = G4Material::GetMaterialTable() ;
393
394 G4double cutinrange ;
395
396
397 G4cout << "give cuts in range" << G4endl ;
398
399 G4cout << "cut for GAMMA in mm =" ;
400 G4cin >> cutinrange ;
401 theGamma->SetCuts(cutinrange) ;
402
403 GammaKineticEnergyCuts = theGamma->GetCutsInEnergy() ;
404
405 G4cout << "cut for ELECTRON in mm =" ;
406 G4cin >> cutinrange ;
407 theElectron->SetCuts(cutinrange) ;
408
409 ElectronKineticEnergyCuts = theElectron->GetCutsInEnergy() ;
410
411 G4cout << "cut for POSITRON in mm =" ;
412 G4cin >> cutinrange ;
413 thePositron->SetCuts(cutinrange) ;
414
415 PositronKineticEnergyCuts = thePositron->GetCutsInEnergy() ;
416
417//****************************************************
418// setcut for the selected particle (if it is not e-/e+)
419 if((theParticle != theElectron) && (theParticle != thePositron))
420 {
421 G4cout << "cut for the selected particle in mm =" ;
422 G4cin >> cutinrange ;
423 theParticle->SetCuts(cutinrange) ;
424
425 ParticleKineticEnergyCuts = theParticle->GetEnergyCuts() ;
426 }
427
428 G4double energy, momentum, mass;
429 G4ProcessVector* palongget ;
430 G4ProcessVector* palongdo ;
431 G4ProcessVector* ppostget ;
432 G4ProcessVector* ppostdo ;
433
434 mass = theParticle->GetPDGMass() ;
435 energy = 1.*GeV + mass ;
436 momentum=std::sqrt(energy*energy-mass*mass) ;
437 G4ParticleMomentum theMomentum(momentum,0.,0.);
438 G4double pModule = theMomentum.mag();
439 G4DynamicParticle aParticle(theParticle,energy,theMomentum);
440 aParticle.SetKineticEnergy(energy-mass);
441
442 outFile << " " << G4endl;
443 outFile << " M S C test **********************************************" << G4endl ;
444 outFile << " " << G4endl;
445 palongget = aParticle.GetDefinition()->GetProcessManager()
446 ->GetAlongStepProcessVector(0);
447 ppostget = aParticle.GetDefinition()->GetProcessManager()
448 ->GetPostStepProcessVector(0);
449 palongdo = aParticle.GetDefinition()->GetProcessManager()
450 ->GetAlongStepProcessVector(1);
451 ppostdo = aParticle.GetDefinition()->GetProcessManager()
452 ->GetPostStepProcessVector(1);
453
454//---------------------------------- Physics --------------------------------
455
456 G4int itry=1, Ntry=1, Nstart, ir;
457 G4double r ;
458
459//**************************************************************************
460 const G4int Nbin=100 ;
461 G4double TkinMeV[Nbin] =
462 {0.0002, 0.0005, 0.001,0.0015,0.002,0.003,0.004,0.005,0.006,0.008,
463 0.01,0.015,0.02,0.03,0.04,0.05,0.06,0.08,
464 0.1,0.15,0.2,0.3,0.4,0.5,0.6,0.8,
465 1.,1.5,2.,3.,4.,5.,6.,8.,
466 10.,15.,20.,30.,40.,50.,60.,80.,
467 100.,150.,200.,300.,400.,500.,600.,800.,
468 1.0e3,1.5e3,2.0e3,3.0e3,4.0e3,5.0e3,6.0e3,8.0e3,
469 1.0e4,1.5e4,2.0e4,3.0e4,4.0e4,5.0e4,6.0e4,8.0e4,
470 1.0e5,1.5e5,2.0e5,3.0e5,4.0e5,5.0e5,6.0e5,8.0e5,
471 1.0e6,1.5e6,2.0e6,3.0e6,4.0e6,5.0e6,6.0e6,8.0e6,
472 1.0e7,1.5e7,2.0e7,3.0e7,4.0e7,5.0e7,6.0e7,8.0e7,
473 1.0e8,1.5e8,2.0e8,3.0e8,4.0e8,5.0e8,6.0e8,8.0e8,
474 1.0e9, 5.0e9} ;
475
476 G4int J=-1 ;
477
478 G4double lambda,trueStep,geomStep,stepLimit,
479 previousStepSize,currentMinimumStep ;
480 G4ParticleChange* aParticleChange ;
481
482 // material loop
483 for (J=0; J<theMaterialTable->length(); J++) {
484
485 apttoMaterial = (*theMaterialTable)[ J ] ;
486 MaterialName = apttoMaterial->GetName() ;
487 G4cout << G4endl << "material=" << MaterialName << G4endl << G4endl;
488
489//---------- Volume definition ---------------------
490
491 G4VPhysicalVolume* myVolume ;
492
493 myVolume = BuildVolume(apttoMaterial) ;
494
495//--------- track and Step definition (for this test ONLY!)------------
496 G4ThreeVector aPosition(0.,0.,0.);
497 const G4ThreeVector aDirection(0.,0.,1.) ;
498 const G4ThreeVector transl(0.,0.,0.) ;
499
500 G4double aTime = 0. ;
501
502 G4Track* tracke = new G4Track(&aParticle,aTime,aPosition) ;
503 G4Track& trackele = (*tracke) ;
504 G4GRSVolume* touche = new G4GRSVolume(myVolume,NULL,transl);
505 (*tracke).SetTouchable(touche);
506
507 (*tracke).SetMomentumDirection(aDirection) ;
508
509 G4Step* Step = new G4Step() ;
510 //******************************************
511 tracke->SetStep(Step);
512 (*Step).InitializeStep(tracke) ;
513
514 const G4Step& Step = (*Step) ;
515
516 G4double CurrentSafety=1000.;
517 G4double& currentSafety = CurrentSafety ;
518
519 G4StepPoint* aPoint = new G4StepPoint();
520 (*aPoint).SetPosition(aPosition) ;
521 G4double safety = 10000.*cm ;
522 (*aPoint).SetSafety(safety) ;
523
524 (*Step).SetPostStepPoint(aPoint) ;
525
526 G4int ib,ev ;
527 G4double xc,yc,zc,later ;
528//**************************************************************************
529
530 //---------------------------------------------------------------------
531 G4double distr[100]=
532 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
533 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
534 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
535 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
536 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
537 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
538 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
539 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
540 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
541 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
542 G4double distrx[100]=
543 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
544 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
545 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
546 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
547 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
548 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
549 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
550 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
551 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
552 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
553
554 G4int nev[100]=
555 {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
556 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
557 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
558 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
559 G4int nevx[100]=
560 {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
561 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
562 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
563 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
564
565 G4double thetamean,thetamean2,dthetamean,errdistr ;
566
567 outFile << G4endl ;
568 outFile << " " << MaterialName << " PostStepDoIt (scattering) test " << G4endl ;
569 outFile << " +++++++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
570 outFile << G4endl ;
571
572 G4double wg[100] ;
573
574 G4double over,overx ;
575 G4double fwg, costheta,theta ;
576 G4int ibin,iw ;
577 G4double TMeV ;
578 G4double dthetadeg,dtheta ;
579 G4int flagdeg ;
580
581//*************************************************************************************
582// example: how to use the GetRange() fuctions............
583
584 for (G4int i=0; i<Nbin; i++) {
585 G4double TMeV= TkinMeV[i]*MeV;
586
587 G4double rangestart;
588 if( theParticle == theElectron )
589 rangestart = theElectronIonisation.OldGetRange(theElectron,TMeV,apttoMaterial) ;
590 if( theParticle == thePositron )
591 rangestart = thePositronIonisation.OldGetRange(thePositron,TMeV,apttoMaterial) ;
592 if( theParticle == theMuonPlus )
593 rangestart = theMuonPlusIonisation.OldGetRange(theMuonPlus,TMeV,apttoMaterial) ;
594 if( theParticle == theMuonMinus )
595 rangestart = theMuonMinusIonisation.OldGetRange(theMuonMinus,TMeV,apttoMaterial) ;
596 if( theParticle == theProton )
597 rangestart = theProtonIonisation.OldGetRange(theProton,TMeV,apttoMaterial) ;
598 if( theParticle == theAntiProton )
599 rangestart = theAntiProtonIonisation.OldGetRange(theAntiProton,TMeV,apttoMaterial) ;
600 if( theParticle == thePionPlus )
601 rangestart = thePionPlusIonisation.OldGetRange(thePionPlus,TMeV,apttoMaterial) ;
602 if( theParticle == thePionMinus )
603 rangestart = thePionMinusIonisation.OldGetRange(thePionMinus,TMeV,apttoMaterial) ;
604 if( theParticle == theKaonPlus )
605 rangestart = theKaonPlusIonisation.OldGetRange(theKaonPlus,TMeV,apttoMaterial) ;
606 if( theParticle == theKaonMinus )
607 rangestart = theKaonMinusIonisation.OldGetRange(theKaonMinus,TMeV,apttoMaterial) ;
608 G4cout << " energy (MeV) = " << TMeV << G4endl;
609 G4cout << " range (OldGetRange) = " << rangestart << G4endl ;
610
611 G4double range2= G4EnergyLossTables::GetRange(theParticle, TMeV, apttoMaterial);
612
613 G4cout << " range (G4EnergyLossTables::GetRange) = " << range2 << G4endl ;
614 }
615
616 } // of material loop
617
618 } // of largest loop
619
620 return EXIT_SUCCESS;
621}
Note: See TracBrowser for help on using the repository browser.