source: trunk/source/processes/electromagnetic/xrays/test/testG4ForwardXrayTR.cc@ 1211

Last change on this file since 1211 was 1199, checked in by garnier, 16 years ago

nvx fichiers dans CVS
File size: 37.3 KB
RevLine 
[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: testG4ForwardXrayTR.cc,v 1.7 2006/06/29 19:56:33 gunter Exp $
28// GEANT4 tag $Name: geant4-09-03-cand-01 $
29//
30
31
32
33#include "G4ios.hh"
34#include <fstream>
35#include <iomanip>
36#include "g4templates.hh"
37#include "globals.hh"
38#include "G4Timer.hh"
39
40#include "G4PhysicsLogVector.hh"
41#include "G4PhysicsFreeVector.hh"
42#include "G4PhysicsVector.hh"
43
44#include "G4eEnergyLoss.hh"
45#include "G4eIonisation.hh"
46#include "G4hEnergyLoss.hh"
47#include "G4PAIenergyLoss.hh"
48#include "G4hIonisation.hh"
49#include "G4PAIonisation.hh"
50#include "G4ForwardXrayTR.hh"
51#include "G4TransitionRadiation.hh"
52
53#include "G4DynamicParticle.hh"
54#include "G4Element.hh"
55#include "G4Material.hh"
56#include "G4PVPlacement.hh"
57#include "G4LogicalVolume.hh"
58#include "G4GRSVolume.hh"
59#include "G4Box.hh"
60#include "G4ProcessManager.hh"
61#include "G4Step.hh"
62#include "G4StepPoint.hh"
63#include "G4Track.hh"
64#include "G4GPILSelection.hh"
65
66#include "G4Electron.hh"
67#include "G4Positron.hh"
68#include "G4Proton.hh"
69#include "G4AntiProton.hh"
70#include "G4PionPlus.hh"
71#include "G4PionMinus.hh"
72#include "G4KaonPlus.hh"
73#include "G4KaonMinus.hh"
74
75// It tests the G4PAIenergyLoss,G4PAIonisation and G4ForwardXrayTR processes
76// created by L.Urban on 06/06/97 and modified by V.Grichine on 17.12.97
77//
78// Modifications:
79// 23-09-97, geometry adapted for the touchable
80// 14.11.97 PAIonisation class is responsible for <dE/dx> calculation
81// 16.12.97 G4ForwardXrayTR class for generation of X-ray TR
82//
83
84G4VPhysicalVolume* BuildVolume(G4Material* matworld)
85// it builds a simple box filled with material matword .......
86{
87 G4Box *myWorldBox= new G4Box ("WBox",10000.*cm,10000.*cm,10000.*cm);
88
89 G4LogicalVolume *myWorldLog = new G4LogicalVolume(myWorldBox,matworld,
90 "WLog",0,0,0) ;
91
92 G4PVPlacement *myWorldPhys = new G4PVPlacement(0,G4ThreeVector(),
93 "WPhys",
94 myWorldLog,
95 0,false,0) ;
96
97
98 return myWorldPhys ;
99
100}
101
102
103int main()
104{
105 //-------- set output format-------
106 G4cout.setf( std::ios::scientific, std::ios::floatfield );
107 //---write results to the file hloss -----
108 std::ofstream outFile("XrayTR.cc", std::ios::out ) ;
109 outFile.setf( std::ios::scientific, std::ios::floatfield );
110
111 //--------- Material definition ---------
112
113 G4double a, z, ez, density ,temperature,pressure;
114 G4State state ;
115 G4String name, symbol;
116 G4int nel;
117
118 G4Element* elH = new G4Element ("Hydrogen", "H", 1. , 1.01*g/mole);
119
120 a = 6.01*g/mole;
121 G4Element* elC = new G4Element(name="Carbon", symbol="C", ez=6., a);
122
123 a = 14.01*g/mole;
124 G4Element* elN = new G4Element(name="Nitrogen", symbol="N", ez=7., a);
125
126 a = 16.00*g/mole;
127 G4Element* elO = new G4Element(name="Oxigen", symbol="O", ez=8., a);
128
129 G4Material* C2H2 = new G4Material ("Polypropelene" , 0.91*g/cm3, 2);
130 C2H2->AddElement(elC,2);
131 C2H2->AddElement(elH,2);
132
133
134 density = 1.29e-03*g/cm3;
135 state = kStateGas ;
136 temperature = 273.*kelvin ;
137 pressure = 1.*atmosphere ;
138 G4Material* Air = new G4Material(name="Air", density, nel=2 ,
139 state ,temperature , pressure ) ;
140 Air->AddElement(elN, .7);
141 Air->AddElement(elO, .3);
142
143
144 /* ***************************************************************
145
146 density = 5.85e-3*g/cm3 ;
147 a = 131.29*g/mole ;
148 G4Material* Xe = new G4Material(name="Xenon", ez=54., a, density);
149
150 density = 1.782e-3*g/cm3 ;
151 a = 39.948*g/mole ;
152 G4Material* Ar = new G4Material(name="Argon", ez=18., a, density);
153
154 a = 9.012*g/mole;
155 density = 1.848*g/cm3;
156 G4Material* Be = new G4Material(name="Beryllium", z=4. , a, density);
157
158 a = 26.98*g/mole;
159 density = 2.7*g/cm3;
160 G4Material* Al = new G4Material(name="Aluminium", z=13., a, density);
161
162 a = 28.09*g/mole;
163 density = 2.33*g/cm3;
164 G4Material* Si = new G4Material(name="Silicon", z=14., a, density);
165
166
167 G4Material* H2O = new G4Material ("Water" , 1.*g/cm3, 2);
168 H2O->AddElement(elH,2);
169 H2O->AddElement(elO,1);
170
171 a = 55.85*g/mole;
172 density = 7.87*g/cm3;
173 G4Material* Fe = new G4Material(name="Iron", z=26., a, density);
174
175 a = 196.97*g/mole;
176 density = 19.32*g/cm3;
177 G4Material* Au = new G4Material(name="Gold", z=79., a, density);
178
179 a = 207.19*g/mole;
180 density = 11.35*g/cm3;
181 G4Material* Pb = new G4Material(name="Lead", z=82., a, density);
182
183 a = 0. ;
184 density = 0. ;
185 G4Material* Vac= new G4Material(name="Vacuum",z=0., a, density,kVacuum);
186
187 ******************************************************************** */
188
189 const G4MaterialTable* theMaterialTable ;
190 G4Material* apttoMaterial ;
191 G4String MaterialName ;
192 G4Timer theTimer ;
193
194
195
196
197//--------- Particle definition ---------
198
199 G4ParticleDefinition* theGamma = G4Gamma::GammaDefinition();
200 G4ParticleDefinition* theElectron = G4Electron::ElectronDefinition();
201 G4ParticleDefinition* thePositron = G4Positron::PositronDefinition();
202 G4ParticleDefinition* theProton = G4Proton::ProtonDefinition();
203 G4ParticleDefinition* theAntiProton = G4AntiProton::AntiProtonDefinition();
204 G4ParticleDefinition* thePionPlus = G4PionPlus::PionPlusDefinition();
205 G4ParticleDefinition* thePionMinus = G4PionMinus::PionMinusDefinition();
206 G4ParticleDefinition* theKaonPlus = G4KaonPlus::KaonPlusDefinition();
207 G4ParticleDefinition* theKaonMinus = G4KaonMinus::KaonMinusDefinition();
208
209 G4double* GammaKineticEnergyCuts ;
210 G4double* ElectronKineticEnergyCuts ;
211 G4double* PositronKineticEnergyCuts ;
212 G4double* ParticleKineticEnergyCuts ;
213
214 theMaterialTable = G4Material::GetMaterialTable() ;
215
216 G4double cutinrange,CutInRangeele,CutInRangepos ;
217
218 G4eIonisation theElectronIonisation , thePositronIonisation ;
219 G4ProcessManager* theElectronProcessManager = theElectron->GetProcessManager();
220 theElectronProcessManager->AddProcess(&theElectronIonisation,-1,0,0) ;
221 G4ProcessManager* thePositronProcessManager = thePositron->GetProcessManager();
222 thePositronProcessManager->AddProcess(&thePositronIonisation,-1,0,0) ;
223
224 G4ParticleWithCuts* theParticle ;
225
226 G4double energy, momentum, mass;
227 G4ProcessVector* palongget ;
228 G4ProcessVector* palongdo ;
229 G4ProcessVector* ppostget ;
230 G4ProcessVector* ppostdo ;
231
232 G4String confirm ;
233
234 G4cout << " Do you want the proton as particle (yes/no)? " << std::flush;
235 //G4cin >> confirm ;
236 confirm = "yes" ;
237 if(confirm == "yes")
238 {
239 mass=theProton->GetPDGMass();
240 theParticle = theProton;
241 }
242 else
243 {
244 G4cout << " Do you want the antiproton as particle (yes/no)? " << std::flush;
245 G4cin >> confirm ;
246 if(confirm == "yes")
247 {
248 mass=theAntiProton->GetPDGMass();
249 theParticle = theAntiProton;
250 }
251 else
252 {
253 G4cout << " Do you want the pi+ as particle (yes/no)? " << std::flush;
254 G4cin >> confirm ;
255 if(confirm == "yes")
256 {
257 mass=thePionPlus->GetPDGMass();
258 theParticle = thePionPlus;
259 }
260 else
261 {
262 G4cout << " Do you want the pi- as particle (yes/no)? " << std::flush;
263 G4cin >> confirm ;
264 if(confirm == "yes")
265 {
266 mass=thePionMinus->GetPDGMass();
267 theParticle = thePionMinus;
268 }
269 else
270 {
271 G4cout << " Do you want the K+ as particle (yes/no)? " << std::flush;
272 G4cin >> confirm ;
273 if(confirm == "yes")
274 {
275 mass=theKaonPlus->GetPDGMass();
276 theParticle = theKaonPlus;
277 }
278 else
279 {
280 G4cout << " Do you want the K- as particle (yes/no)? " << std::flush;
281 G4cin >> confirm ;
282 if(confirm == "yes")
283 {
284 mass=theKaonMinus->GetPDGMass();
285 theParticle = theKaonMinus;
286 }
287 else
288 {
289 G4cout << " There is no other particle in the test." << G4endl;
290 return EXIT_FAILURE;
291 }
292 }
293 }
294 }
295 }
296 }
297
298
299 energy = 1000.0*GeV + mass ; // was 1.0*GeV now 1.0*TeV
300 momentum=std::sqrt(energy*energy-mass*mass) ;
301
302 G4ParticleMomentum theMomentum(momentum,0.,0.);
303
304 G4double pModule = theMomentum.mag();
305
306 G4DynamicParticle aParticle(theParticle,energy,theMomentum);
307
308 aParticle.SetKineticEnergy(energy-mass);
309
310
311 // G4hIonisation theParticleIonisation ;
312
313 G4PAIonisation theParticleIonisation ;
314
315 G4ProcessManager* theParticleProcessManager = theParticle->GetProcessManager();
316 theParticleProcessManager->AddProcess(&theParticleIonisation,-1,0,0) ;
317
318 G4ForceCondition cond ;
319 G4ForceCondition* condition = &cond ;
320
321 G4double currentSafety ;
322 G4double& refsafety=currentSafety;
323
324
325 G4cout << "cut for GAMMA in mm =" ;
326 // G4cin >> cutinrange ;
327 cutinrange = 0.1 ;
328 theGamma->SetCuts(cutinrange) ;
329 G4cout << "gamma,cut in range(mm)=" << theGamma->GetCuts() << G4endl ;
330 outFile << " ---------------------------------------" << G4endl ;
331 outFile << " gamma,cut in range(mm)=" << theGamma->GetCuts() << G4endl ;
332
333 GammaKineticEnergyCuts = theGamma->GetCutsInEnergy() ;
334 for (G4int icut=0; icut<theMaterialTable->length(); icut++)
335 {
336 G4cout << "material index=" << icut << " kin.energy cut(MeV)=" <<
337 GammaKineticEnergyCuts[icut] << G4endl ;
338 outFile << " material index=" << icut << " kin.energy cut(MeV)=" <<
339 GammaKineticEnergyCuts[icut] << G4endl ;
340 }
341
342
343 G4cout << "cut for ELECTRON in mm =" ;
344 // G4cin >> cutinrange ;
345 cutinrange = 1.0 ;
346 CutInRangeele = cutinrange ;
347 theElectron->SetCuts(cutinrange) ;
348 G4cout << "electron,cut in range(mm)=" << theElectron->GetCuts() << G4endl ;
349 outFile << " ---------------------------------------" << G4endl ;
350 outFile << " electron,cut in range(mm)=" << theElectron->GetCuts() << G4endl ;
351
352 ElectronKineticEnergyCuts = theElectron->GetCutsInEnergy() ;
353 for ( icut=0; icut<theMaterialTable->length(); icut++)
354 {
355 G4cout << "material index=" << icut << " kin.energy cut(MeV)=" <<
356 ElectronKineticEnergyCuts[icut] << G4endl ;
357 outFile << " material index=" << icut << " kin.energy cut(MeV)=" <<
358 ElectronKineticEnergyCuts[icut] << G4endl ;
359 }
360
361 G4cout << "cut for POSITRON in mm =" ;
362 // G4cin >> cutinrange ;
363 cutinrange = 1.0 ;
364 CutInRangepos = cutinrange ;
365 thePositron->SetCuts(cutinrange) ;
366 G4cout << "positron,cut in range(mm)=" << thePositron->GetCuts() << G4endl ;
367 outFile << " ---------------------------------------" << G4endl ;
368 outFile << " positron,cut in range(mm)=" << thePositron->GetCuts() << G4endl ;
369
370 PositronKineticEnergyCuts = thePositron->GetCutsInEnergy() ;
371 for ( icut=0; icut<theMaterialTable->length(); icut++)
372 {
373 G4cout << "material index=" << icut << " kin.energy cut(MeV)=" <<
374 PositronKineticEnergyCuts[icut] << G4endl ;
375 outFile << " material index=" << icut << " kin.energy cut(MeV)=" <<
376 PositronKineticEnergyCuts[icut] << G4endl ;
377 }
378
379 G4cout << "cut for hadrons in mm =" ;
380 // G4cin >> cutinrange ;
381 cutinrange = 1.0 ;
382 theParticle->SetCuts(cutinrange) ;
383 G4cout << "after particle setcuts " << G4endl;
384
385
386 G4cout << "cut in range(mm)=" << theParticle->GetLengthCuts() << G4endl ;
387 outFile << " ---------------------------------------" << G4endl ;
388 outFile << " cut in range(mm)=" << theParticle->GetLengthCuts() << G4endl ;
389
390 ParticleKineticEnergyCuts = theParticle->GetEnergyCuts() ;
391
392 for ( icut=0; icut<theMaterialTable->length(); icut++)
393 {
394 G4cout << "material index=" << icut << " kin.energy cut(MeV)=" <<
395 ParticleKineticEnergyCuts[icut] << G4endl ;
396 outFile << " material index=" << icut << " kin.energy cut(MeV)=" <<
397 ParticleKineticEnergyCuts[icut] << G4endl ;
398 }
399
400 G4cout << " ------ ----- " << G4endl ;
401 outFile << " " << G4endl;
402
403 // Creation of X-ray transition radiation object
404 G4cout<<"Just before TR object creation"<<G4endl ;
405
406 G4ForwardXrayTR theXrayTR ;
407
408 G4cout<<"Just after TR object creation"<<G4endl ;
409 theParticleProcessManager->AddProcess(&theXrayTR,-1,0,0) ;
410
411
412
413
414 /* *********************************************************************
415 // Output to file of fPAItransferBank data
416
417 G4bool isOutRange ;
418 G4PhysicsLogVector*
419 aLogVector = new G4PhysicsLogVector(G4PAIonisation::GetMinKineticEnergy(),
420 G4PAIonisation::GetMaxKineticEnergy(),
421 G4PAIonisation::GetBinNumber() ) ;
422
423 for(G4int materialIndex=0 ;
424 materialIndex < theMaterialTable->length() ; materialIndex++)
425 {
426 apttoMaterial = (*theMaterialTable)[materialIndex] ;
427 MaterialName = apttoMaterial->GetName() ;
428 outFile<<"PAI transfer data for the material = "
429 <<MaterialName<<G4endl<<G4endl ;
430
431 outFile<<"Particle Tkin"<<"\t\t"
432 <<"Transfer energy, keV"<<"\t"
433 <<"Primary ionisation, 1/mm"<<G4endl<<G4endl ;
434
435 // G4cout<<"Read of TotBin ? TotBin = "<<G4PAIonisation::GetBinNumber()<<G4endl;
436
437 for(G4int iTkin=48;iTkin<G4PAIonisation::GetBinNumber();iTkin++)
438 {
439
440 outFile<<aLogVector->GetLowEdgeEnergy(iTkin)<<G4endl<<G4endl ;
441
442 G4PhysicsVector* trVec = (*G4PAIenergyLoss::GetPAItransferBank())
443 (materialIndex*G4PAIonisation::GetBinNumber() + iTkin) ;
444 G4cout<<materialIndex*G4PAIonisation::GetBinNumber() + iTkin
445 <<"\t\t"<<trVec->GetVectorLength()<<G4endl ;
446 for(G4int iTr=0;iTr<trVec->GetVectorLength();iTr++)
447 {
448 outFile<<"\t\t"<<trVec->GetLowEdgeEnergy(iTr)*1000.
449 <<"\t\t"<<(*trVec)(iTr)<<G4endl ;
450 }
451
452 }
453 outFile<<G4endl ;
454 }
455
456
457 G4cout<<"Start dE/dx distribution"<<G4endl ;
458 G4int iLoss, iStat, iStatMax = 10000;
459 G4double energyLoss[200], delta ;
460 G4int spectrum[200] ;
461 const G4DynamicParticle* particlePtr = &aParticle ;
462
463 for(iLoss=0;iLoss<200;iLoss++)
464 {
465 energyLoss[iLoss] = 0.5*iLoss*keV ;
466 spectrum[iLoss] = 0 ;
467 }
468 for(iStat=0;iStat<iStatMax;iStat++)
469 {
470 delta = G4PAIenergyLoss::GetLossWithFluct(10.0*mm,particlePtr,Xe) ;
471 for(iLoss=0;iLoss<200;iLoss++)
472 {
473 if(delta <= energyLoss[iLoss]) break ;
474 }
475 spectrum[iLoss-1]++ ;
476 }
477 G4double meanLoss = 0.0 ;
478 outFile<<"Energy loss, keV"<<"\t\t"<<"Distribution"<<G4endl<<G4endl ;
479 for(iLoss=0;iLoss<200;iLoss++)
480 {
481 outFile<<energyLoss[iLoss]*1000.0<<"\t\t"<<spectrum[iLoss]<<G4endl ;
482 meanLoss +=energyLoss[iLoss]*spectrum[iLoss] ;
483 }
484 G4cout<<"Mean loss over spectrum = "<<meanLoss*1000./iStatMax<<" keV"<<G4endl ;
485
486
487 **************************************************************** */
488
489 // Output of TR information
490
491 G4int iMat, jMat, iPlace, iTkin, iTR ;
492 G4PhysicsLogVector* vectorTkin = new
493 G4PhysicsLogVector(G4ForwardXrayTR::GetMinProtonTkin(),
494 G4ForwardXrayTR::GetMaxProtonTkin(),
495 G4ForwardXrayTR::GetTotBin() ) ;
496
497 for(iMat=0;iMat<theMaterialTable->length();iMat++)
498 {
499 for(jMat=0;jMat<theMaterialTable->length();jMat++)
500 {
501 if(jMat == iMat) continue ;
502
503 if(jMat < iMat)
504 {
505 outFile<<(*theMaterialTable)[iMat]->GetName()<<" -> "<<"\t\t"
506 <<(*theMaterialTable)[jMat]->GetName()<<G4endl<<G4endl ;
507 outFile<<"Lorentz Factor"<<"\t\t"
508 <<"energy TR number"<<"\t\t"
509 <<"theta TR number"<<"\t\t"<<G4endl<<G4endl ;
510 for(iTkin=0;iTkin<G4ForwardXrayTR::GetTotBin();iTkin++)
511 {
512 iPlace =(iMat*(theMaterialTable->length()-1)+jMat)*
513 G4ForwardXrayTR::GetTotBin()+iTkin ;
514
515 outFile<<1+(vectorTkin->GetLowEdgeEnergy(iTkin)/proton_mass_c2)<<"\t\t"
516 <<(*(*G4ForwardXrayTR::GetEnergyDistrTable())
517 (iPlace))(0)<<"\t\t"
518 <<(*(*G4ForwardXrayTR::GetAngleDistrTable())
519 (iPlace))(0)<<G4endl ;
520 G4cout<<1+(vectorTkin->GetLowEdgeEnergy(iTkin)/proton_mass_c2)<<"\t\t"
521 <<(*(*G4ForwardXrayTR::GetEnergyDistrTable())
522 (iPlace))(0)<<"\t\t"
523 <<(*(*G4ForwardXrayTR::GetAngleDistrTable())
524 (iPlace))(0)<<G4endl ;
525 }
526 }
527 else
528 {
529 outFile<<(*theMaterialTable)[iMat]->GetName()<<" -> "<<"\t\t"
530 <<(*theMaterialTable)[jMat]->GetName()<<G4endl<<G4endl ;
531 outFile<<"Lorentz Factor"<<"\t\t"
532 <<"energy TR number"<<"\t\t"
533 <<"theta TR number"<<"\t\t"<<G4endl<<G4endl ;
534 for(iTkin=0;iTkin<G4ForwardXrayTR::GetTotBin();iTkin++)
535 {
536 iPlace =(iMat*(theMaterialTable->length()-1)+jMat-1)*
537 G4ForwardXrayTR::GetTotBin()+iTkin ;
538
539 outFile<<1+(vectorTkin->GetLowEdgeEnergy(iTkin)/proton_mass_c2)<<"\t\t"
540 <<(*(*G4ForwardXrayTR::GetEnergyDistrTable())
541 (iPlace))(0)<<"\t\t"
542 <<(*(*G4ForwardXrayTR::GetAngleDistrTable())
543 (iPlace))(0)<<G4endl ;
544 G4cout<<1+(vectorTkin->GetLowEdgeEnergy(iTkin)/proton_mass_c2)<<"\t\t"
545 <<(*(*G4ForwardXrayTR::GetEnergyDistrTable())
546 (iPlace))(0)<<"\t\t"
547 <<(*(*G4ForwardXrayTR::GetAngleDistrTable())
548 (iPlace))(0)<<G4endl ;
549 }
550 }
551 outFile<<G4endl ;
552 G4cout<<G4endl ;
553 }
554 }
555
556 for(iMat=0;iMat<theMaterialTable->length();iMat++)
557 {
558 for(jMat=0;jMat<theMaterialTable->length();jMat++)
559 {
560 if(jMat == iMat) continue ;
561
562 if(jMat < iMat)
563 {
564 outFile<<(*theMaterialTable)[iMat]->GetName()<<" -> "<<"\t\t"
565 <<(*theMaterialTable)[jMat]->GetName()<<G4endl<<G4endl ;
566 for(iTkin=0;iTkin<G4ForwardXrayTR::GetTotBin();iTkin++)
567 {
568 if(iTkin == 9 || iTkin == 19 || iTkin == 29
569 || iTkin == 39 || iTkin == 49 )
570 {
571 outFile<<"iTkin = "<<iTkin<<"\t\t"
572 <<"Lorentz factor = "
573 <<1+(vectorTkin->GetLowEdgeEnergy(iTkin)/proton_mass_c2)
574 <<G4endl<<G4endl ;
575 outFile<<"TR Energy"<<"\t\t"
576 <<" > energy TR number"<<"\t\t"
577 <<"TR ThetaSq"<<"\t\t"
578 <<" > theta TR number"<<"\t\t"<<G4endl<<G4endl ;
579 iPlace =(iMat*(theMaterialTable->length()-1)+jMat)*
580 G4ForwardXrayTR::GetTotBin()+iTkin ;
581
582 for(iTR=0;iTR<G4ForwardXrayTR::GetBinTR();iTR++)
583 {
584 outFile<<(*G4ForwardXrayTR::GetEnergyDistrTable())
585 (iPlace)->GetLowEdgeEnergy(iTR)<<"\t\t"
586 <<(*(*G4ForwardXrayTR::GetEnergyDistrTable())
587 (iPlace))(iTR)<<"\t\t"
588 <<(*G4ForwardXrayTR::GetAngleDistrTable())
589 (iPlace)->GetLowEdgeEnergy(iTR)<<"\t\t"
590 <<(*(*G4ForwardXrayTR::GetAngleDistrTable())
591 (iPlace))(iTR)<<G4endl ;
592 }
593 }
594 }
595 }
596 else
597 {
598 outFile<<(*theMaterialTable)[iMat]->GetName()<<" -> "<<"\t\t"
599 <<(*theMaterialTable)[jMat]->GetName()<<G4endl<<G4endl ;
600 for(iTkin=0;iTkin<G4ForwardXrayTR::GetTotBin();iTkin++)
601 {
602 if(iTkin == 9 || iTkin == 19 || iTkin == 29
603 || iTkin == 39 || iTkin == 49 )
604 {
605 outFile<<"iTkin = "<<iTkin<<"\t\t"
606 <<"Lorentz factor = "
607 <<1+(vectorTkin->GetLowEdgeEnergy(iTkin)/proton_mass_c2)
608 <<G4endl<<G4endl ;
609 outFile<<"TR Energy"<<"\t\t"
610 <<" > energy TR number"<<"\t\t"
611 <<"TR ThetaSq"<<"\t\t"
612 <<" > theta TR number"<<"\t\t"<<G4endl<<G4endl ;
613 iPlace =(iMat*(theMaterialTable->length()-1)+jMat-1)*
614 G4ForwardXrayTR::GetTotBin()+iTkin ;
615
616 for(iTR=0;iTR<G4ForwardXrayTR::GetBinTR();iTR++)
617 {
618 outFile<<(*G4ForwardXrayTR::GetEnergyDistrTable())
619 (iPlace)->GetLowEdgeEnergy(iTR)<<"\t\t"
620 <<(*(*G4ForwardXrayTR::GetEnergyDistrTable())
621 (iPlace))(iTR)<<"\t\t"
622 <<(*G4ForwardXrayTR::GetAngleDistrTable())
623 (iPlace)->GetLowEdgeEnergy(iTR)<<"\t\t"
624 <<(*(*G4ForwardXrayTR::GetAngleDistrTable())
625 (iPlace))(iTR)<<G4endl ;
626 }
627 }
628 }
629 }
630 outFile<<G4endl ;
631 G4cout<<G4endl ;
632 }
633 }
634
635 G4cout<<"Start TR energy distribution"<<G4endl ;
636 G4int iLossTR, iStatTR, iStatMaxTR = 1000;
637 G4double energyTR[200], deltaTR ;
638 G4int spectrumTR[200] ;
639 // const G4DynamicParticle* particlePtr = &aParticle ;
640
641 for(iLossTR=0;iLossTR<200;iLossTR++)
642 {
643 energyTR[iLossTR] = 0.5*iLossTR*keV ;
644 spectrumTR[iLossTR] = 0 ;
645 }
646 for(iStatTR=0;iStatTR<iStatMaxTR;iStatTR++)
647 {
648 deltaTR = theXrayTR.GetEnergyTR(0,1,39) ;
649 for(iLossTR=0;iLossTR<200-1;iLossTR++)
650 {
651 if(deltaTR <= energyTR[iLossTR]) break ;
652 }
653 spectrumTR[iLossTR]++ ;
654 }
655 G4double meanLossTR = 0.0 ;
656 outFile<<"TR energy, keV"<<"\t\t"<<"Distribution"<<G4endl<<G4endl ;
657 for(iLossTR=0;iLossTR<200;iLossTR++)
658 {
659 outFile<<energyTR[iLossTR]*1000.0<<"\t\t"<<spectrumTR[iLossTR]<<G4endl ;
660 meanLossTR +=energyTR[iLossTR]*spectrumTR[iLossTR] ;
661 }
662 G4cout<<"Mean TR energy over spectrum = "
663 <<meanLossTR*1000./iStatMaxTR<<" keV"<<G4endl ;
664
665
666 G4int message ;
667 G4cout<<"Continue ?(enter int>0) Break ? (enter int<=0)"<<G4endl ;
668 G4cin>>message ;
669 if(message<= 0) // break program at this point
670 {
671 return 1 ;
672 }
673
674
675
676 outFile << " ionisation test **************************************" << G4endl ;
677 outFile << " " << G4endl;
678 outFile << " particle = " <<
679 aParticle.GetDefinition()->GetParticleName() << G4endl ;
680 outFile << G4endl;
681
682 palongget = aParticle.GetDefinition()->GetProcessManager()
683 ->GetAlongStepProcessVector(typeGPIL);
684 ppostget = aParticle.GetDefinition()->GetProcessManager()
685 ->GetPostStepProcessVector(typeGPIL);
686 palongdo = aParticle.GetDefinition()->GetProcessManager()
687 ->GetAlongStepProcessVector(typeDoIt);
688 ppostdo = aParticle.GetDefinition()->GetProcessManager()
689 ->GetPostStepProcessVector(typeDoIt);
690
691//---------------------------------- Physics --------------------------------
692
693 G4int itry=1, Ntry=1, Nstart, ir;
694 G4double r ;
695
696//**************************************************************************
697 const G4int Nbin=97 ;
698 G4double TkinMeV[Nbin] =
699 {0.001,0.0015,0.002,0.003,0.004,0.005,0.006,0.008,
700 0.01,0.015,0.02,0.03,0.04,0.05,0.06,0.08,
701 0.1,0.15,0.2,0.3,0.4,0.5,0.6,0.8,
702 1.,1.5,2.,3.,4.,5.,6.,8.,
703 10.,15.,20.,30.,40.,50.,60.,80.,
704 100.,150.,200.,300.,400.,500.,600.,800.,
705 1.0e3,1.5e3,2.0e3,3.0e3,4.0e3,5.0e3,6.0e3,8.0e3,
706 1.0e4,1.5e4,2.0e4,3.0e4,4.0e4,5.0e4,6.0e4,8.0e4,
707 1.0e5,1.5e5,2.0e5,3.0e5,4.0e5,5.0e5,6.0e5,8.0e5,
708 1.0e6,1.5e6,2.0e6,3.0e6,4.0e6,5.0e6,6.0e6,8.0e6,
709 1.0e7,1.5e7,2.0e7,3.0e7,4.0e7,5.0e7,6.0e7,8.0e7,
710 1.0e8,1.5e8,2.0e8,3.0e8,4.0e8,5.0e8,6.0e8,8.0e8,
711 1.0e9} ;
712
713 G4int J=-1 ;
714
715 G4double lambda,trueStep,geomStep,stepLimit,
716 previousStepSize,currentMinimumStep ;
717 G4ParticleChange* aParticleChange ;
718
719 G4double T,dEdx,range ;
720
721 NEXTMATERIAL: ;
722 J = J+1 ;
723 if ( J >= theMaterialTable->length() )
724 { G4cout << "that was the last material in the table --> STOP" << G4endl;
725 return EXIT_FAILURE ; }
726
727 apttoMaterial = (*theMaterialTable)[ J ] ;
728 MaterialName = apttoMaterial->GetName() ;
729 G4cout << "material=" << MaterialName << G4endl ;
730 G4cout << "Do you want the Energyloss test 1. for this material?" << G4endl ;
731 G4cout << "type a positive number if the answer is YES" << G4endl ;
732 G4cout << "type a negative number if the answer is NO " << G4endl ;
733 G4int icont ;
734 G4cin >> icont ;
735 if ( icont < 0 )
736 goto NEXTMATERIAL ;
737
738//---------- Volume definition ---------------------
739
740 G4VPhysicalVolume* myVolume ;
741
742 myVolume = BuildVolume(apttoMaterial) ;
743
744//--------- track and Step definition (for this test ONLY!)------------
745 G4ThreeVector aPosition(0.,0.,0.);
746 const G4ThreeVector aDirection(0.,0.,1.) ;
747 G4double aTime = 0. ;
748
749 G4Track* tracke = new G4Track(&aParticle,aTime,aPosition) ;
750 G4Track& trackele = (*tracke) ;
751 //(*tracke).SetVolume(myVolume) ;
752 G4GRSVolume* touche = new G4GRSVolume(myVolume, NULL, aPosition);
753 (*tracke).SetTouchable(touche);
754 (*tracke).SetMomentumDirection(aDirection) ;
755
756
757 G4Step* Step = new G4Step() ;
758 G4Step& Step = (*Step) ;
759
760 G4StepPoint* aPoint = new G4StepPoint();
761 (*aPoint).SetPosition(aPosition) ;
762 G4double safety = 10000.*cm ;
763 (*aPoint).SetSafety(safety) ;
764
765 (*Step).SetPostStepPoint(aPoint) ;
766
767//**************************************************************************
768
769 G4cout << G4endl;
770 G4cout <<" " << MaterialName << " Energyloss test 1." << G4endl ;
771 G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
772 G4cout << G4endl ;
773 G4cout << "kin.en.(MeV) dE/dx(MeV/mm) range(mm) Step(mm)" << G4endl ;
774 G4cout << G4endl ;
775
776 outFile << G4endl;
777 outFile <<" " << MaterialName << " Energyloss test 1." << G4endl ;
778 outFile << " +++++++++++++++++++++++++++++++++++++++++" << G4endl ;
779 outFile << G4endl ;
780 outFile << "kin.en.(MeV) dE/dx(MeV/mm) range(mm) Step(mm)" << G4endl ;
781 outFile << G4endl ;
782
783
784 G4GPILSelection* selection ;
785
786 for ( G4int i=0 ; i<Nbin ; i++)
787 {
788 trueStep = cutinrange ;
789 previousStepSize = cutinrange ;
790 currentMinimumStep = trueStep ;
791 (*tracke).SetKineticEnergy(TkinMeV[i]) ;
792 stepLimit = (*palongget)(0)->AlongStepGetPhysicalInteractionLength(
793 trackele,
794 previousStepSize,
795 currentMinimumStep,
796 refsafety,
797 selection) ;
798
799 dEdx = theParticleIonisation.GetdEdx() ;
800 range = theParticleIonisation.GetRangeNow() ;
801 T = TkinMeV[i] ;
802
803 G4cout <<" " << T << " " << dEdx/mm << " " ;
804 G4cout << range/mm << " " << stepLimit/mm << G4endl ;
805
806 outFile <<" " << T << " " << dEdx/mm << " " ;
807 outFile << range/mm << " " << stepLimit/mm << G4endl ;
808
809 }
810
811 G4cout << G4endl;
812 outFile << G4endl;
813
814 ENERGYLOSS2: ;
815
816 G4cout << "material=" << MaterialName << G4endl ;
817 G4cout << "Do you want the Energyloss test 2. for this material?" << G4endl ;
818 G4cout << "type a positive number if the answer is YES" << G4endl ;
819 G4cout << "type a negative number if the answer is NO " << G4endl ;
820 G4cin >> icont ;
821 if ( icont < 0 )
822 goto ENERGYLOSS3 ;
823
824 G4double TMeV,stepmm,stepmx,meanloss,lossnow ;
825
826
827 G4cout << "give an energy value in MeV " ;
828 G4cin >> TMeV ;
829
830 trueStep = cutinrange ;
831 previousStepSize = cutinrange ;
832 currentMinimumStep = trueStep ;
833 (*tracke).SetKineticEnergy(TMeV) ;
834 stepmx = (*palongget)(0)->AlongStepGetPhysicalInteractionLength(
835 trackele,
836 previousStepSize,
837 currentMinimumStep,
838 refsafety,
839 selection);
840
841 G4cout << " give a steplength in mm , the max. meaningful Step is " << stepmx << " mm" <<G4endl;
842 G4cout << "Step:" ;
843 G4cin >> stepmm ;
844
845 (*Step).SetTrack(tracke) ;
846 (*Step).SetStepLength(stepmm);
847
848
849 aParticleChange = (G4ParticleChange*)((*palongdo)(0)->
850 AlongStepDoIt(trackele,Step));
851 meanloss = theParticleIonisation.GetMeanLoss() ;
852 lossnow = TMeV-(*aParticleChange).GetEnergyChange();
853
854 G4cout << G4endl;
855 G4cout <<" " << MaterialName << " Energyloss test 2." << G4endl ;
856 G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
857 G4cout << G4endl ;
858 G4cout << "kin.en.(MeV) Step(mm) meanloss(MeV) act.loss(MeV)" << G4endl ;
859 G4cout << TMeV << " " << stepmm << " " << meanloss << " " << lossnow << G4endl ;
860 G4cout << " status change:" << (*aParticleChange).GetStatusChange() << G4endl ;
861 G4cout << G4endl ;
862
863 outFile << G4endl;
864 outFile <<" " << MaterialName << " Energyloss test 2." << G4endl ;
865 outFile << " +++++++++++++++++++++++++++++++++++++++++" << G4endl ;
866 outFile << G4endl ;
867 outFile << "kin.en.(MeV) Step(mm) meanloss(MeV) act.loss(MeV)" << G4endl ;
868 outFile << TMeV << " " << stepmm << " " << meanloss << " " << lossnow << G4endl ;
869 outFile << " status change:" << (*aParticleChange).GetStatusChange() << G4endl ;
870 outFile << G4endl ;
871
872 goto ENERGYLOSS2 ;
873
874 ENERGYLOSS3: ;
875
876 G4cout << "material=" << MaterialName << G4endl ;
877 G4cout << "Do you want the Energyloss test 3. for this material?" << G4endl ;
878 G4cout << "type a positive number if the answer is YES" << G4endl ;
879 G4cout << "type a negative number if the answer is NO " << G4endl ;
880 G4cin >> icont ;
881 if ( icont < 0 )
882 goto DELTARAY1 ;
883
884 G4cout << "give an energy value in MeV " ;
885 G4cin >> TMeV ;
886
887 trueStep = cutinrange ;
888 previousStepSize = cutinrange ;
889 currentMinimumStep = trueStep ;
890 (*tracke).SetKineticEnergy(TMeV) ;
891 stepmx = (*palongget)(0)->AlongStepGetPhysicalInteractionLength(
892 trackele,
893 previousStepSize,
894 currentMinimumStep,
895 refsafety,
896 selection);
897
898 G4cout << " give a steplength in mm , the max. meaningful Step is " << stepmx << " mm" <<G4endl;
899 G4cout << "Step:" ;
900 G4cin >> stepmm ;
901
902 (*Step).SetTrack(tracke) ;
903 (*Step).SetStepLength(stepmm);
904
905
906 G4cout << " give number of events you want " ;
907 G4int nbev,ibev ;
908 G4cin >> nbev ;
909
910 meanloss=0.;
911 theTimer.Start();
912
913 for ( ibev=0; ibev<nbev; ibev++)
914 {
915 aParticleChange =(G4ParticleChange*)((*palongdo)(0)->
916 AlongStepDoIt(trackele,Step));
917 lossnow = TMeV-(*aParticleChange).GetEnergyChange();
918
919 meanloss += lossnow ;
920 }
921
922 theTimer.Stop();
923 meanloss /= nbev ;
924 G4cout << G4endl;
925 G4cout <<" " << MaterialName << " Energyloss test 3." << G4endl ;
926 G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
927 G4cout << G4endl ;
928 G4cout << "kin.en.(MeV) Step(mm) meanloss(MeV) time/event(sec) " << G4endl ;
929 G4cout << TMeV << " " << stepmm << " " << meanloss << " " <<
930 theTimer.GetUserElapsed()/nbev << G4endl ;
931 G4cout << G4endl ;
932
933 outFile << G4endl;
934 outFile <<" " << MaterialName << " Energyloss test 3." << G4endl ;
935 outFile << " +++++++++++++++++++++++++++++++++++++++++" << G4endl ;
936 outFile << G4endl ;
937 outFile << "kin.en.(MeV) Step(mm) meanloss(MeV) time/event(sec) " << G4endl ;
938 outFile << TMeV << " " << stepmm << " " << meanloss << " " <<
939 theTimer.GetUserElapsed()/nbev << G4endl ;
940 outFile << G4endl ;
941
942 goto ENERGYLOSS3 ;
943
944 DELTARAY1: ;
945 G4cout << "material=" << MaterialName << G4endl ;
946 G4cout << "Do you want the delta ray test 1. for this material?" << G4endl ;
947 G4cout << "type a positive number if the answer is YES" << G4endl ;
948 G4cout << "type a negative number if the answer is NO " << G4endl ;
949 G4cin >> icont ;
950 if ( icont < 0 )
951 goto DELTARAY2 ;
952
953
954 G4cout << G4endl;
955 G4cout <<" " << MaterialName << " delta ray test 1." << G4endl ;
956 G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
957 G4cout << G4endl ;
958 G4cout << "kin.en.(MeV) mean free path(mm)" << G4endl ;
959 G4cout << G4endl ;
960
961 outFile << G4endl;
962 outFile <<" " << MaterialName << " delta ray test 1." << G4endl ;
963 outFile << " +++++++++++++++++++++++++++++++++++++++++" << G4endl ;
964 outFile << G4endl ;
965 outFile << "kin.en.(MeV) mean free path(mm)" << G4endl ;
966 outFile << G4endl ;
967
968 for ( i=0 ; i<Nbin ; i++)
969 {
970
971 previousStepSize = cutinrange ;
972 (*tracke).SetKineticEnergy(TkinMeV[i]) ;
973 stepLimit = theParticleIonisation.GetMeanFreePath(
974 trackele,
975 previousStepSize,
976 condition) ;
977
978 T = TkinMeV[i] ;
979
980 G4cout <<" " << T << " " << stepLimit/mm << G4endl ;
981
982 outFile <<" " << T << " " << stepLimit/mm << G4endl ;
983
984 }
985
986 G4cout << G4endl;
987 outFile << G4endl;
988
989//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
990 DELTARAY2: ;
991
992 G4cout << "material=" << MaterialName << G4endl ;
993 G4cout << "Do you want the deltaray test 2. for this material?" << G4endl ;
994 G4cout << "type a positive number if the answer is YES" << G4endl ;
995 G4cout << "type a negative number if the answer is NO " << G4endl ;
996 G4cin >> icont ;
997
998 G4double newenergy,dx,dy,dz,Tdelta,ddx,ddy,ddz ;
999 G4int nd ;
1000 const G4ThreeVector* momdir ;
1001 G4ParticleMomentum ddir ;
1002
1003 if ( icont < 0 )
1004 goto BREMS1 ;
1005
1006 G4cout << "give an energy value in MeV " ;
1007 G4cin >> TMeV ;
1008
1009 stepmm = 1. ;
1010
1011 (*Step).SetTrack(tracke) ;
1012 (*Step).SetStepLength(stepmm);
1013
1014
1015 (*tracke).SetKineticEnergy(TMeV) ;
1016 aParticleChange = (G4ParticleChange*)((*ppostdo)(0)->
1017 PostStepDoIt(trackele,Step));
1018
1019 newenergy=(*aParticleChange).GetEnergyChange() ;
1020 momdir=(*aParticleChange).GetMomentumChange();
1021 dx = (*momdir).x();
1022 dy = (*momdir).y();
1023 dz = (*momdir).z();
1024 nd=aParticleChange->GetNumberOfSecondaries();
1025
1026 if(nd>0)
1027 {
1028 Tdelta=aParticleChange->GetSecondary(0)->GetKineticEnergy();
1029 ddir=aParticleChange->GetSecondary(0)->
1030 GetMomentumDirection();
1031 ddx = (ddir).x();
1032 ddy = (ddir).y();
1033 ddz = (ddir).z();
1034 }
1035
1036 G4cout << G4endl;
1037 G4cout <<" " << MaterialName << " delta ray test 2." << G4endl ;
1038 G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
1039 G4cout << G4endl ;
1040 G4cout << "T=" << TMeV << " newT=" << newenergy << " (MeV)" << G4endl ;
1041 G4cout << " status change:" << (*aParticleChange).GetStatusChange() << G4endl ;
1042 if(nd>0)
1043 G4cout << "Tdelta=" << Tdelta << G4endl ;
1044 G4cout << "new direction:" << dx << " " << dy << " " << dz << G4endl;
1045 if(nd>0)
1046 G4cout << "delta direction:" << ddx << " " << ddy << " " << ddz << G4endl ;
1047 G4cout << G4endl ;
1048
1049 outFile << G4endl;
1050 outFile <<" " << MaterialName << " delta ray test 2." << G4endl ;
1051 outFile << " +++++++++++++++++++++++++++++++++++++++++" << G4endl ;
1052 outFile << G4endl ;
1053 outFile << "T=" << TMeV << " newT=" << newenergy << " (MeV)" << G4endl;
1054 outFile << " status change:" << (*aParticleChange).GetStatusChange() << G4endl ;
1055 if(nd>0)
1056 outFile << "Tdelta=" << Tdelta << G4endl ;
1057 outFile << "new direction:" << dx << " " << dy << " " << dz << G4endl;
1058 if(nd>0)
1059 outFile << "delta direction:" << ddx << " " << ddy << " " << ddz << G4endl ;
1060 outFile << G4endl ;
1061
1062 (*aParticleChange).Clear();
1063
1064 goto DELTARAY2 ;
1065
1066 BREMS1: ;
1067
1068 if( J < theMaterialTable->length()-1 )
1069 goto NEXTMATERIAL ;
1070
1071 return EXIT_SUCCESS;
1072}
Note: See TracBrowser for help on using the repository browser.