source: trunk/source/processes/electromagnetic/lowenergy/test/G4LowEnergyPolarizedRayleighTest.cc@ 1350

Last change on this file since 1350 was 1350, checked in by garnier, 15 years ago

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