source: trunk/source/processes/electromagnetic/lowenergy/test/G4DNATest.cc@ 1201

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

nvx fichiers dans CVS
File size: 27.6 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// $Id: G4DNATest.cc,v 1.14 2006/06/29 19:43:57 gunter Exp $
27// GEANT4 tag $Name: geant4-09-03-cand-01 $
28
29#include "globals.hh"
30#include "G4ios.hh"
31#include <fstream>
32#include <iomanip>
33#include <memory>
34#include <cstdlib>
35
36#include "G4ParticleDefinition.hh"
37#include "G4ParticleTypes.hh"
38#include "G4ParticleTable.hh"
39#include "G4Material.hh"
40#include "G4MaterialTable.hh"
41#include "G4ProcessManager.hh"
42
43#include "G4EnergyLossTables.hh"
44#include "G4VParticleChange.hh"
45#include "G4ParticleChange.hh"
46#include "G4DynamicParticle.hh"
47#include "G4ForceCondition.hh"
48
49#include "G4LowEnergyBremsstrahlung.hh"
50#include "G4LowEnergyIonisation.hh"
51#include "G4eIonisation.hh"
52#include "G4MultipleScattering.hh"
53#include "G4eIonisation.hh"
54#include "G4eBremsstrahlung.hh"
55#include "G4eplusAnnihilation.hh"
56
57#include "G4ComptonScattering.hh"
58#include "G4PhotoElectricEffect.hh"
59
60#include "G4RunManager.hh"
61
62#include "G4Electron.hh"
63#include "G4Positron.hh"
64#include "G4Gamma.hh"
65
66#include "G4GRSVolume.hh"
67#include "G4Box.hh"
68#include "G4PVPlacement.hh"
69#include "G4Step.hh"
70#include "G4ProductionCutsTable.hh"
71#include "G4MaterialCutsCouple.hh"
72
73#include "G4UnitsTable.hh"
74
75#include "AIDA/IManagedObject.h"
76#include "AIDA/IAnalysisFactory.h"
77#include "AIDA/ITreeFactory.h"
78#include "AIDA/ITree.h"
79#include "AIDA/IHistogramFactory.h"
80#include "AIDA/IHistogram1D.h"
81#include "AIDA/IHistogram2D.h"
82#include "AIDA/IHistogram3D.h"
83#include "AIDA/ITupleFactory.h"
84#include "AIDA/ITuple.h"
85
86// DNA
87#include "G4DNAGenericIonsManager.hh"
88#include "G4DNAElectronElasticBrenner.hh"
89#include "G4DNAElectronElasticEmfietzoglou.hh"
90#include "G4DNAProtonExcitation.hh"
91#include "G4DNAHeliumExcitation.hh"
92#include "G4DNAAlphaPlusExcitation.hh"
93#include "G4DNAAlphaPlusPlusExcitation.hh"
94#include "G4DNAElectronExcitation.hh"
95#include "G4DNAProtonRuddIonization.hh"
96#include "G4DNAProtonChargeDecrease.hh"
97#include "G4DNAHydrogenChargeIncrease.hh"
98#include "G4DNAHydrogenRuddIonization.hh"
99#include "G4DNAProtonBornExcitation.hh"
100#include "G4DNAElectronBornExcitation.hh"
101#include "G4DNAAlphaPlusChargeDecrease.hh"
102#include "G4DNAAlphaPlusChargeIncrease.hh"
103#include "G4DNAAlphaPlusPlusChargeDecrease.hh"
104#include "G4DNAHeliumChargeIncrease.hh"
105
106//! \brief Options structure
107struct Options
108{
109 //! \brief Mean free path test
110 bool meanFreePathTest;
111 //! \brief Post step do it test
112 bool postStepDoItTest;
113 //! \brief Post step do it test
114 bool randomEnergy;
115 //! \brief Output file name
116 const char *outputFileName;
117 //! \brief Material name
118 const char *material;
119 //! \brief Process name
120 const char *process;
121 //! \brief Particle name
122 const char *particle;
123 //! \brief Minimum energy
124 G4double minEnergy;
125 //! \brief Maximum energy
126 G4double maxEnergy;
127 //! \brief Number of energy step
128 G4int nEnergySteps;
129 //! \brief Number of interactions
130 G4int nIterations;
131};
132
133//! \brief Default output file name
134struct Options defaultOptions = { false, false, false, "G4DNATest.hbook", "Water", "G4DNAElectronElasticBrenner", "electron", 7.5*eV, 200*eV, 300, 1 };
135
136//! \brief Creates some materials
137void CreateMaterials(void)
138{
139 G4Element * H = new G4Element ("Hydrogen", "H", 1., 1.01*g/mole);
140 G4Element * O = new G4Element ("Oxygen", "O", 8., 16.00*g/mole);
141 G4Element * C = new G4Element ("Carbon", "C", 6., 12.00*g/mole);
142 G4Element * Cs = new G4Element ("Cesium", "Cs", 55., 132.905*g/mole);
143 G4Element * I = new G4Element ("Iodine", "I", 53., 126.9044*g/mole);
144
145 G4Material * Si = new G4Material("Silicon", 14., 28.055*g/mole, 2.33*g/cm3);
146 G4Material * Fe = new G4Material("Iron", 26., 55.85*g/mole, 7.87*g/cm3);
147 G4Material * Cu = new G4Material("Copper", 29., 63.55*g/mole, 8.96*g/cm3);
148 G4Material * W = new G4Material("Tungsten", 74., 183.85*g/mole, 19.30*g/cm3);
149 G4Material * Pb = new G4Material("Lead", 82., 207.19*g/mole, 11.35*g/cm3);
150 G4Material * U = new G4Material("Uranium", 92., 238.03*g/mole, 18.95*g/cm3);
151 G4Material * maO = new G4Material("Oxygen", 8., 16.00*g/mole, 1.1*g/cm3);
152 G4Material * water = new G4Material ("Water", 1.*g/cm3, 2);
153 water->AddElement(H, 2);
154 water->AddElement(O, 1);
155
156 G4Material* ethane = new G4Material ("Ethane", 0.4241*g/cm3, 2);
157 ethane->AddElement(H, 6);
158 ethane->AddElement(C, 2);
159
160 G4Material* csI = new G4Material ("CsI", 4.53*g/cm3, 2);
161 csI->AddElement(Cs, 1);
162 csI->AddElement(I, 1);
163
164 // This is needed to suppress some warnings. These lines can be deleted;
165 Si->GetTemperature();
166 Fe->GetTemperature();
167 Cu->GetTemperature();
168 W->GetTemperature();
169 Pb->GetTemperature();
170 U->GetTemperature();
171 maO->GetTemperature();
172 water->GetTemperature();
173 ethane->GetTemperature();
174 csI->GetTemperature();
175}
176
177//! \brief Process the options arguments
178//! \param argc Number of arguments
179//! \param argv Pointer to the arguments
180//! \param options Structure to fill-in
181void processOptions(int argc, char ** argv, struct Options * options)
182{
183 options->meanFreePathTest = defaultOptions.meanFreePathTest;
184 options->postStepDoItTest = defaultOptions.meanFreePathTest;
185 options->randomEnergy = defaultOptions.randomEnergy;
186 options->outputFileName = defaultOptions.outputFileName;
187 options->material = defaultOptions.material;
188 options->process = defaultOptions.process;
189 options->particle = defaultOptions.particle;
190 options->minEnergy = defaultOptions.minEnergy;
191 options->maxEnergy = defaultOptions.maxEnergy;
192 options->nEnergySteps = defaultOptions.nEnergySteps;
193 options->nIterations = defaultOptions.nIterations;
194
195 int i(1);
196
197 while (i<argc)
198 {
199 if (argv[i][0]=='-' && argv[i][2]==0)
200 {
201 switch(argv[i][1])
202 {
203 case 'h':
204 case '?':
205 G4cout << argv[0] << " [-h|-?] [-a] [-b] [-r] [-o <file name>] [-m <material name>] [-p <process name>] [-P <particle name>] [-e <min energy in eV>] [-E <max energy in eV>] [-s <energy steps>] [-n <iterations>] " << G4endl
206 << G4endl
207 << "-h|-? Shows this help" << G4endl
208 << "-a Enables mean free path test" << G4endl
209 << "-b Enables post step do it test" << G4endl
210 << "-r Energy is choosen at random within the range" << G4endl
211 << "-o <arg> Set the output file name (default: \"" << defaultOptions.outputFileName << "\")" << G4endl
212 << "-m <arg> Set the material (default: \"" << defaultOptions.material << "\")" << G4endl
213 << "-p <arg> Set the process (default: \"" << defaultOptions.process << "\")" << G4endl
214 << "-P <arg> Set the incoming particle (default: \"" << defaultOptions.particle << "\")" << G4endl
215 << "-e <arg> Set the low energy range in eV (default: " << defaultOptions.minEnergy/eV << " eV)" << G4endl
216 << "-E <arg> Set the high energy range in eV (default: " << defaultOptions.maxEnergy/eV << " eV)" << G4endl
217 << "-s <arg> Set the energy range step (default: " << defaultOptions.nEnergySteps << ")"<< G4endl
218 << "-n <arg> Set the number of iterations for the post step do it (default: " << defaultOptions.nIterations << ")" << G4endl;
219 exit(0);
220 break;
221
222 case 'a':
223 options->meanFreePathTest=true;
224 break;
225
226 case 'b':
227 options->postStepDoItTest=true;
228 break;
229
230 case 'r':
231 options->randomEnergy=true;
232 break;
233
234 case 'o':
235 i++;
236 if (i<argc)
237 {
238 options->outputFileName = argv[i];
239 break;
240 }
241
242 case 'm':
243 i++;
244 if (i<argc)
245 {
246 options->material = argv[i];
247 break;
248 }
249
250 case 'p':
251 i++;
252 if (i<argc)
253 {
254 options->process = argv[i];
255 break;
256 }
257
258 case 'P':
259 i++;
260 if (i<argc)
261 {
262 options->particle = argv[i];
263 break;
264 }
265
266 case 'e':
267 i++;
268 if (i<argc)
269 {
270 options->minEnergy = std::atof(argv[i])*eV;
271 if (options->minEnergy <= 0.)
272 {
273 G4cout << argv[0] << ": Energy must be > 0." << G4endl;
274 exit(-1);
275 }
276
277 break;
278 }
279
280 case 'E':
281 i++;
282 if (i<argc)
283 {
284 options->maxEnergy = std::atof(argv[i])*eV;
285 if (options->maxEnergy <= 0.)
286 {
287 G4cout << argv[0] << ": Energy must be > 0." << G4endl;
288 exit(-1);
289 }
290
291 break;
292 }
293
294 case 's':
295 i++;
296 if (i<argc)
297 {
298 options->nEnergySteps = atoi(argv[i]);
299 if (options->nEnergySteps <= 1)
300 {
301 G4cout << argv[0] << ": Expected at least two steps." << G4endl;
302 exit(-1);
303 }
304
305 break;
306 }
307
308 G4cout << argv[0] << ": Expected one more parameter in " << argv[i] << " option. Use -h option for help." << G4endl;
309 exit(-1);
310
311 case 'n':
312 i++;
313 if (i<argc)
314 {
315 options->nIterations = atoi(argv[i]);
316 if (options->nIterations <= 0)
317 {
318 G4cout << argv[0] << ": Expected at least one iteration." << G4endl;
319 exit(-1);
320 }
321
322 break;
323 }
324
325 G4cout << argv[0] << ": Expected one more parameter in " << argv[i] << " option. Use -h option for help." << G4endl;
326 exit(-1);
327
328 default:
329 G4cout << argv[0] << ": Unknown " << argv[i] << " option. Use -h option for help." << G4endl;
330 exit(-1);
331 }
332 }
333 else
334 {
335 G4cout << argv[0] << ": Bad arguments. Use -h option for help." << G4endl;
336 exit(-1);
337 }
338
339 i++;
340 }
341
342 if (options->minEnergy >= options->maxEnergy)
343 {
344 G4cout << argv[0] << ": Mininum energy is higher than maximum energy" << G4endl;
345 exit(-1);
346 }
347
348 G4cout << "Mean free path test: ";
349 if (options->meanFreePathTest)
350 G4cout << "On";
351 else
352 G4cout << "Off";
353 G4cout << G4endl << "Post step do it test: ";
354 if (options->postStepDoItTest)
355 G4cout << "On";
356 else
357 G4cout << "Off";
358 G4cout << G4endl << "Random energy generation: ";
359 if (options->randomEnergy)
360 G4cout << "On";
361 else
362 G4cout << "Off";
363 G4cout << G4endl << "Output file: " << options->outputFileName << G4endl;
364 G4cout << "Material: " << options->material << G4endl;
365 G4cout << "Process: " << options->process << G4endl;
366 G4cout << "Min energy: " << options->minEnergy/eV << " eV" << G4endl;
367 G4cout << "Max energy: " << options->maxEnergy/eV << " eV" << G4endl;
368 G4cout << "N energy steps: " << options->nEnergySteps << G4endl;
369 G4cout << "N iterations: " << options->nIterations << G4endl;
370}
371
372//! \brief Return the selected material
373//! \param options Options for the material choice
374//! \return The material
375G4Material * GetSelectedMaterial(const struct Options & options)
376{
377 const G4MaterialTable* theMaterialTable=G4Material::GetMaterialTable();
378
379 G4int i(G4Material::GetNumberOfMaterials());
380
381 while (i>0)
382 {
383 i--;
384
385 if ((*theMaterialTable)[i]->GetName()==options.material)
386 return (*theMaterialTable)[i];
387 }
388
389 i=G4Material::GetNumberOfMaterials();
390
391 G4cout << "Available materials are: " << G4endl;
392 while (i>0)
393 {
394 i--;
395 G4cout << (*theMaterialTable)[i]->GetName();
396
397 if (i>0)
398 G4cout << ", ";
399 }
400
401 G4cout << G4endl;
402
403 exit(-2);
404 return 0;
405}
406
407//! \brief Creates the geometry
408//! \param options Options for the material choice
409//! \return The world volume
410G4PVPlacement * CreateGeometry(const struct Options & options)
411{
412 G4Box* theFrame = new G4Box ("Frame", 1*mm, 1*mm, 1*mm);
413
414 G4LogicalVolume* logicalFrame = new G4LogicalVolume(theFrame, GetSelectedMaterial(options), "LFrame", 0, 0, 0);
415
416 G4PVPlacement * placement = new G4PVPlacement(0, G4ThreeVector(), "PFrame", logicalFrame, 0, false, 0);
417
418 G4cout << "[OK] Geometry built" << G4endl;
419 return placement;
420}
421
422//! \brief Get process from options
423//! \param options Options for the process choice
424//! \return The choosen process
425G4VLowEnergyTestableDiscreteProcess * GetSelectedProcess(const struct Options & options)
426{
427 static G4VLowEnergyTestableDiscreteProcess ** processes=0;
428 if (!processes)
429 {
430 G4DNAGenericIonsManager * genericIonsManager;
431 genericIonsManager=G4DNAGenericIonsManager::Instance();
432
433 processes=new G4VLowEnergyTestableDiscreteProcess * [18];
434 processes[0]=new G4DNAElectronElasticBrenner;
435 processes[1]=new G4DNAElectronElasticEmfietzoglou;
436 processes[2]=new G4DNAProtonExcitation;
437 processes[3]=new G4DNAHeliumExcitation;
438 processes[4]=new G4DNAAlphaPlusExcitation;
439 processes[5]=new G4DNAAlphaPlusPlusExcitation;
440 processes[6]=new G4DNAElectronExcitation;
441 processes[7]=new G4DNAProtonRuddIonization;
442 processes[8]=new G4DNAProtonChargeDecrease;
443 processes[9]=new G4DNAHydrogenChargeIncrease;
444 processes[10]=new G4DNAHydrogenRuddIonization;
445 processes[11]=new G4DNAProtonBornExcitation;
446 processes[12]=new G4DNAElectronBornExcitation;
447 processes[13]=new G4DNAAlphaPlusChargeDecrease;
448 processes[14]=new G4DNAAlphaPlusChargeIncrease;
449 processes[15]=new G4DNAAlphaPlusPlusChargeDecrease;
450 processes[16]=new G4DNAHeliumChargeIncrease;
451 processes[17]=0;
452 }
453
454 unsigned long i(0);
455 while (processes[i])
456 {
457 if (processes[i]->GetProcessName()==options.process)
458 return processes[i];
459
460 i++;
461 }
462
463 G4cout << "Available processes are: " << G4endl;
464 i=0;
465 while (processes[i])
466 {
467 G4cout << processes[i]->GetProcessName();
468 i++;
469
470 if (processes[i])
471 G4cout << ", ";
472 }
473
474 G4cout << G4endl;
475
476 exit(-2);
477 return 0;
478}
479
480G4ParticleDefinition * const * ParticleList()
481{
482 static G4ParticleDefinition * * particles=0;
483
484 if (!particles)
485 {
486 G4DNAGenericIonsManager * genericIonsManager;
487 genericIonsManager=G4DNAGenericIonsManager::Instance();
488
489 particles=new G4ParticleDefinition * [10];
490 particles[0]=G4Electron::Electron();
491 particles[1]=G4Positron::Positron();
492 particles[2]=G4Gamma::Gamma();
493 particles[3]=G4AntiProton::AntiProton();
494 particles[4]=G4Proton::Proton();
495 particles[5]=genericIonsManager->GetIon("alpha++");
496 particles[6]=genericIonsManager->GetIon("alpha+");
497 particles[7]=genericIonsManager->GetIon("helium");
498 particles[8]=genericIonsManager->GetIon("hydrogen");
499 particles[9]=0;
500 }
501
502 return particles;
503}
504
505G4int GetParticleIndex(G4ParticleDefinition * const particle)
506{
507 G4int i(0);
508 G4ParticleDefinition * const * particles(ParticleList());
509
510 while (*particles)
511 {
512 if ((*particles)==particle)
513 return i;
514
515 i++;
516 particles++;
517 }
518
519 return -1;
520}
521
522G4ParticleDefinition * GetSelectedParticle(const struct Options & options)
523{
524 G4ParticleDefinition * const * particles(ParticleList());
525 while (*particles)
526 {
527 if ((*particles)->GetParticleName()==options.particle)
528 return (*particles);
529
530 particles++;
531 }
532
533 G4cout << "Available particles are: " << G4endl;
534 particles=ParticleList();
535
536 while (*particles)
537 {
538 G4cout << (*particles)->GetParticleName();
539
540 particles++;
541
542 if (*particles)
543 G4cout << ", ";
544 }
545
546 G4cout << G4endl;
547
548 exit(-2);
549 return 0;
550}
551
552//! \brief Setup processes
553//! \param options Options for the process choice
554void SetPhysics(const struct Options & options, G4VPhysicalVolume * world)
555{
556 G4ProductionCutsTable * cutsTable(G4ProductionCutsTable::GetProductionCutsTable());
557 G4ProductionCuts * cuts(cutsTable->GetDefaultProductionCuts());
558
559 G4ParticleDefinition * const * particles(ParticleList());
560 while (*particles)
561 {
562 cuts->SetProductionCut(1*micrometer, *particles);
563 particles++;
564 }
565
566 cutsTable->UpdateCoupleTable(world);
567 G4cout << "[OK] Cuts are defined " << G4endl;
568
569 G4ParticleDefinition * particle(GetSelectedParticle(options));
570 G4VProcess * eProcess=GetSelectedProcess(options);
571 if (! (eProcess->IsApplicable(*particle)))
572 {
573 G4cout<< "Process " << eProcess->GetProcessName() << " is not applicable to " << particle->GetParticleName() << G4endl;
574 exit(0);
575 return;
576 }
577
578 G4cout<< "[OK] Process " << eProcess->GetProcessName() << " is applicable to " << particle->GetParticleName() << G4endl;
579
580 G4ProcessManager * gProcessManager(new G4ProcessManager(particle));
581 particle->SetProcessManager(gProcessManager);
582 gProcessManager->AddDiscreteProcess(eProcess);
583
584 G4cout << "[OK] Processes are defined " << G4endl;
585
586
587 G4cout << "[OK] Building physics tables" << G4endl;
588 eProcess->BuildPhysicsTable(* particle);
589
590 G4cout << "[OK] Physics tables built" << G4endl;
591}
592
593//! \brief Generates the step
594//! \param options Options related to the track generation
595//! \return The generated track
596G4Step * GenerateStep(const struct Options & options)
597{
598 G4ThreeVector momentumDirection;
599
600 momentumDirection.setRThetaPhi(1., std::acos(2*G4UniformRand()-1.), twopi * G4UniformRand());
601
602 G4double lnEnergyMin=std::log(options.minEnergy);
603 G4double lnEnergyMax=std::log(options.maxEnergy);
604 G4DynamicParticle * dynamicElectron(new G4DynamicParticle(GetSelectedParticle(options), momentumDirection, std::exp(lnEnergyMin+(lnEnergyMax-lnEnergyMin)*G4UniformRand())));
605
606 G4Track * aTrack(new G4Track(dynamicElectron, 0., G4ThreeVector(0., 0., 0.)));
607
608 G4Step* aStep(new G4Step());
609 aStep->SetTrack(aTrack);
610 aTrack->SetStep(aStep);
611
612 G4Material * material(GetSelectedMaterial(options));
613 G4ProductionCutsTable * cutsTable(G4ProductionCutsTable::GetProductionCutsTable());
614 const G4MaterialCutsCouple * theCouple(cutsTable->GetMaterialCutsCouple(material, cutsTable->GetDefaultProductionCuts()));
615
616 G4StepPoint * aPoint(new G4StepPoint());
617 aPoint->SetPosition(G4ThreeVector(0., 0., 0.));
618 aPoint->SetMaterial(material);
619 aPoint->SetMaterialCutsCouple(theCouple);
620 aPoint->SetSafety(10000.*cm);
621
622 aStep->SetPreStepPoint(aPoint);
623
624 return aStep;
625}
626
627void ProgressBar(G4int remainingIterations)
628{
629 static time_t startingTime;
630 static time_t nextDumpTime;
631 static G4int startingIteration(0);
632 time_t now;
633
634 if (remainingIterations==0)
635 {
636 startingIteration=0;
637 }
638 else if (startingIteration==0)
639 {
640 startingTime=time(0);
641 nextDumpTime=startingTime+3;
642 startingIteration=remainingIterations;
643 }
644 else
645 {
646 now=time(0);
647 if (now>nextDumpTime)
648 {
649 nextDumpTime=now+10;
650 G4double time;
651 G4double perc;
652
653 time=std::floor(static_cast<G4double>(now-startingTime)/static_cast<G4double>(startingIteration-remainingIterations)*static_cast<G4double>(remainingIterations)+0.5);
654 perc=std::floor(static_cast<G4double>(remainingIterations)/static_cast<G4double>(startingIteration)*200.+.5)/2.;
655
656 G4cout << " " << perc << " % Remaining time: " << time << " s \r";
657 G4cout.flush();
658 }
659 }
660}
661
662//! \brief Test the mean free path table
663//! \param tupleFactory The tuple factory
664//! \param options Options related to the mean free path test
665void MeanFreePathTest(AIDA::ITupleFactory * tupleFactory, const struct Options & options)
666{
667 AIDA::ITuple* iTuple = tupleFactory->create("1", "Mean Free Path Ntuple", "double k, log_k, mfp, log_mfp, cpu_time");
668
669 G4double energy(options.minEnergy);
670 G4double stpEnergy(std::pow(options.maxEnergy/energy, 1./static_cast<G4double>(options.nEnergySteps-1)));
671 G4int step(options.nEnergySteps);
672
673 G4ForceCondition condition;
674 G4VLowEnergyTestableDiscreteProcess * process(GetSelectedProcess(options));
675
676 G4double mfp;
677 clock_t time;
678
679 ProgressBar(0);
680 while (step>0)
681 {
682 G4Step * aStep(GenerateStep(options));
683 G4Track * aTrack(aStep->GetTrack());
684
685 if (!options.randomEnergy)
686 {
687 aTrack->SetKineticEnergy(energy);
688 energy*=stpEnergy;
689 }
690 ProgressBar(step);
691 step--;
692
693 time=clock();
694 mfp=process->DumpMeanFreePath(*aTrack, 1.*mm, &condition)/m;
695 time=clock()-time;
696
697 iTuple->fill(iTuple->findColumn("k"), aTrack->GetKineticEnergy()/eV);
698 iTuple->fill(iTuple->findColumn("log_k"), std::log10(aTrack->GetKineticEnergy()/eV));
699 iTuple->fill(iTuple->findColumn("mfp"), mfp);
700 iTuple->fill(iTuple->findColumn("log_mfp"), std::log10(mfp));
701 iTuple->fill(iTuple->findColumn("cpu_time"), static_cast<G4double>(time)/static_cast<G4double>(CLOCKS_PER_SEC));
702 iTuple->addRow();
703
704 delete aTrack;
705 delete aStep;
706 }
707}
708
709//! \brief Test the post step do it
710//! \param tupleFactory The tuple factory
711//! \param options Options related to the post step do it test
712void PostStepDoItTest(AIDA::ITupleFactory * tupleFactory, const struct Options & options)
713{
714 AIDA::ITuple* iTuple = tupleFactory->create("2", "Post Step Do It Test", "double iteration, step, in_k, log_in_k, in_theta, in_phi, in_pol_theta, in_pol_phi, e_deposit, log_e_deposit, trk_status, out_k, log_out_k, out_theta, out_phi, out_pol_theta, out_pol_phi, n_secondaries, sec0_type, sec0_k, log_sec0_k, sec0_theta, sec0_phi, sec0_pol_theta, sec0_pol_phi, sec1_type, sec1_k, log_sec1_k, sec1_theta, sec1_phi, sec1_pol_theta, sec1_pol_phi, sec2_type, sec2_k, log_sec2_k, sec2_theta, sec2_phi, sec2_pol_theta, sec2_pol_phi, sec3_type, sec3_k, log_sec3_k, sec3_theta, sec3_phi, sec3_pol_theta, sec3_pol_phi, cpu_time");
715
716 G4double energy(options.minEnergy);
717 G4double stpEnergy(std::pow(options.maxEnergy/energy, 1./static_cast<G4double>(options.nEnergySteps-1)));
718 G4int step(options.nEnergySteps);
719
720 G4VLowEnergyTestableDiscreteProcess * process(GetSelectedProcess(options));
721 clock_t time;
722
723 char strBuffer[20];
724
725 ProgressBar(0);
726 while (step>0)
727 {
728 G4Step * aStep(GenerateStep(options));
729 G4Track * aTrack(aStep->GetTrack());
730 const G4DynamicParticle * aParticle(aTrack->GetDynamicParticle());
731 G4ThreeVector vector;
732
733 if (!options.randomEnergy)
734 {
735 aTrack->SetKineticEnergy(energy);
736 energy*=stpEnergy;
737 }
738 ProgressBar(step);
739 step--;
740
741 G4int iteration(options.nIterations);
742
743 while (iteration>0)
744 {
745 iteration--;
746
747 aStep->SetStepLength(1*micrometer);
748
749 iTuple->fill(iTuple->findColumn("iteration"), iteration);
750 iTuple->fill(iTuple->findColumn("step"), aStep->GetStepLength()/m);
751
752 iTuple->fill(iTuple->findColumn("in_k"), aParticle->GetKineticEnergy()/eV);
753 iTuple->fill(iTuple->findColumn("log_in_k"), std::log10(aParticle->GetKineticEnergy()/eV));
754 vector=aParticle->GetMomentumDirection();
755 iTuple->fill(iTuple->findColumn("in_theta"), vector.theta());
756 iTuple->fill(iTuple->findColumn("in_phi"), vector.phi());
757 vector=aParticle->GetPolarization();
758 iTuple->fill(iTuple->findColumn("in_pol_theta"), vector.theta());
759 iTuple->fill(iTuple->findColumn("in_pol_phi"), vector.phi());
760
761 time=clock();
762 G4ParticleChange * particleChange(dynamic_cast<G4ParticleChange *>(process->PostStepDoIt(*aTrack, *aStep)));
763 time=clock()-time;
764
765 aTrack->SetKineticEnergy(particleChange->GetEnergy());
766 aTrack->SetMomentumDirection(*particleChange->GetMomentumDirection());
767 aTrack->SetPolarization(*particleChange->GetPolarization());
768
769 iTuple->fill(iTuple->findColumn("e_deposit"), particleChange->GetLocalEnergyDeposit()/eV);
770 iTuple->fill(iTuple->findColumn("log_e_deposit"), std::log10(particleChange->GetLocalEnergyDeposit()/eV));
771 iTuple->fill(iTuple->findColumn("trk_status"), particleChange->GetTrackStatus());
772
773 iTuple->fill(iTuple->findColumn("out_k"), aParticle->GetKineticEnergy()/eV);
774 iTuple->fill(iTuple->findColumn("log_out_k"), std::log10(aParticle->GetKineticEnergy()/eV));
775 vector=aParticle->GetMomentumDirection();
776 iTuple->fill(iTuple->findColumn("out_theta"), vector.theta());
777 iTuple->fill(iTuple->findColumn("out_phi"), vector.phi());
778 vector=aParticle->GetPolarization();
779 iTuple->fill(iTuple->findColumn("out_pol_theta"), vector.theta());
780 iTuple->fill(iTuple->findColumn("out_pol_phi"), vector.phi());
781
782 G4int n(particleChange->GetNumberOfSecondaries());
783 iTuple->fill(iTuple->findColumn("n_secondaries"), n);
784
785 while (n>0)
786 {
787 n--;
788 G4Track * aSecTrack(particleChange->GetSecondary(n));
789 const G4DynamicParticle * aSecParticle(aSecTrack->GetDynamicParticle());
790
791 sprintf(strBuffer, "sec%d_type", n);
792 iTuple->fill(iTuple->findColumn(strBuffer), GetParticleIndex(aSecParticle->GetDefinition()));
793 sprintf(strBuffer, "sec%d_k", n);
794 iTuple->fill(iTuple->findColumn(strBuffer), aSecParticle->GetKineticEnergy()/eV);
795 sprintf(strBuffer, "log_sec%d_k", n);
796 iTuple->fill(iTuple->findColumn(strBuffer), std::log10(aSecParticle->GetKineticEnergy()/eV));
797 vector=aSecParticle->GetMomentumDirection();
798 sprintf(strBuffer, "sec%d_theta", n);
799 iTuple->fill(iTuple->findColumn(strBuffer), vector.theta());
800 sprintf(strBuffer, "sec%d_phi", n);
801 iTuple->fill(iTuple->findColumn(strBuffer), vector.phi());
802 vector=aSecParticle->GetPolarization();
803 sprintf(strBuffer, "sec%d_pol_theta", n);
804 iTuple->fill(iTuple->findColumn(strBuffer), vector.theta());
805 sprintf(strBuffer, "sec%d_pol_phi", n);
806 iTuple->fill(iTuple->findColumn(strBuffer), vector.phi());
807
808 delete aSecTrack;
809 }
810
811 iTuple->fill(iTuple->findColumn("cpu_time"), static_cast<G4double>(time)/static_cast<G4double>(CLOCKS_PER_SEC));
812 iTuple->addRow();
813
814 particleChange->Clear();
815 }
816
817 delete aTrack;
818 delete aStep;
819 }
820}
821
822//! \brief Main function
823//! \param argc Number of arguments
824//! \param argv Pointer to the arguments
825//! \return The exit value
826int main(int argc, char ** argv)
827{
828 struct Options options;
829 processOptions(argc, argv, &options);
830
831 CreateMaterials();
832
833 GetSelectedProcess(options);
834 GetSelectedMaterial(options);
835
836 G4RunManager* rm = new G4RunManager();
837 rm->GeometryHasBeenModified();
838 G4VPhysicalVolume * world(CreateGeometry(options));
839 rm->DefineWorldVolume(world);
840 G4cout << "[OK] World is defined " << G4endl;
841
842 SetPhysics(options, world);
843
844 if (!(options.meanFreePathTest || options.postStepDoItTest))
845 {
846 G4cout << "[OK] Program completed" << G4endl;
847 return 0;
848 }
849
850 // HBOOK initialization
851 AIDA::IAnalysisFactory * analysisFactory(AIDA_createAnalysisFactory());
852 AIDA::ITreeFactory * treeFactory(analysisFactory->createTreeFactory());
853 AIDA::ITree * tree(treeFactory->create(options.outputFileName, "hbook", false, true));
854 G4cout << "[OK] Tree store: " << tree->storeName() << G4endl;
855
856 AIDA::ITupleFactory * tupleFactory(analysisFactory->createTupleFactory(*tree));
857
858 // Mean free path test
859 if (options.meanFreePathTest)
860 {
861 G4cout << "[OK] Mean free path test started" << G4endl;
862 MeanFreePathTest(tupleFactory, options);
863 G4cout << "[OK] Mean free path test completed" << G4endl;
864 }
865
866 // Post step do it test
867 if (options.postStepDoItTest)
868 {
869 G4cout << "[OK] Post step do it test started" << G4endl;
870 PostStepDoItTest(tupleFactory, options);
871 G4cout << "[OK] Post step do it test completed" << G4endl;
872 }
873
874 G4cout << "[OK] Storing analysis data" << G4endl;
875 tree->commit();
876 tree->close();
877
878 G4cout << "[OK] Deleting analysis data" << G4endl;
879 delete tupleFactory;
880 delete tree;
881 delete treeFactory;
882 delete analysisFactory;
883
884 G4cout << "[OK] Program completed" << G4endl;
885 return 0;
886}
Note: See TracBrowser for help on using the repository browser.