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
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: 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
70G4VPhysicalVolume* 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
89int 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(energy*energy-mass*mass) ;
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.